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ABSTRACT 

The development of gravitational wave observatories (Advanced LIGO/Virgo, Einstein Telescope) is 
proceeding apace, and the direct detection of gravitational waves should be imminent. The last decade 
of observational and theoretical developments in stellar and binary evolution provides an opportunity 
to incorporate major improvements to the predictions from populations synthesis models. We update 
the StarTrack code, and compute the rates for NS-NS, BH-NS, and BH-BH mergers. Among the most 
important revisions in the formation and evolution of double compact binaries are: updated wind mass 
loss rates (allowing for stellar mass black holes up to 80 M©), a realistic treatment of the common 
envelope phase (a process that can affect merger rates by 2-3 orders of magnitude), and a qualitatively 
new neutron star/black hole mass distribution (consistent with the observed "mass gap"). We present 
a parameter study with these major physical updates included, focusing on the most important factors 
that set the double compact object merger rates. A few of our more interesting findings include: (i) 
The binding energy of the envelope plays a pivotal role in determining whether a binary merges within 
a Hubble time, and our results depend sensitively on the model describing this, (ii) Our description of 
natal kicks from supernovae plays an important role, especially for the formation of BH-BH systems. 
(in) The masses of BH-BH systems can be substantially increased in the case of low mctallicities or 
weak winds, (iv) Certain combinations of parameters underpredict the Galactic NS-NS merger rate, 
and can be ruled out. (v) Models incorporating delayed supernovae do not agree with the observed 
NS/BH "mass gap", in accordance with our previous work. And, finally, (vi) We find enhanced rates for 
BH-BH mergers as compared to previous estimates, with an expectation of ^ 100 such mergers per year 
in Advanced LIGO/Virgo detectors (although this rate is sensitive to the details of underlying models, 
such as the natal kick distribution). This is the first in a series of three papers. The second paper will 
study the merger rates of double compact objects as a function of cosmological redshift, star formation 
rate, and metallicity. In the third paper we will present the detection rates for future gravitational wave 
observatories, using up-to-date signal waveforms and sensitivity curves. 
Subject headings: common envelope, supernova, black hole, neutron star 



1. INTRODUCTION 

We investigate the evolution of binary systems that 
leads to the formation of close double compact objects 
(DCOs) merging within a Hubble time: double neutron 
star (NS-NS), black hole-neutron star (BH-NS), and dou- 
ble black hole (BH-BH) systems. Only isolated evolution 
(i.e., in field populations) is considered. Several different 
groups have provided s imilar estimates in the past decade 
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[20031: iNutzman et aU 120041: [Pfahl et all . I2005D , and explo- 
ration of the relevant physical processes and rate estimates 
were obtained with the StarTrack 
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I2010allbl . [2011a: ,0'Shaughnessv et al.l . ,2005a,b,y). In this 
paper we revisit these rate estimates, incorporating recent 
improvements to the input physics within the StarTrack 
code. Similar rate estimates performed for dense pop- 
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ulations in which dynamical interactions between stars 
are importa nt (i.e., globula r cluste r s) have been pr e sented 



elsewher e ([Giiltekin et al, 
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To describe the formation of a double compact object 
one needs to know the history of both progenitors, as well 
as the interactions between them. However, the uncertain- 
ties associated with a number of evolutionary processes in- 
hibit an accurate description. Among the most important 
unknowns are the details of the common envelope (CE) 
phase, the unknown maximum mass of the NS, the physics 
of the supernova explosions that form compact objects, 
and the wind mass loss rates of their progenitors. In this 
paper we approach the problem by calculating and investi- 
gating a suite of population synthesis models. The uncer- 
tainties are assessed by altering the parameters and input 
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physics within defined hmits. This ahows us to "bracket" 
our ignorance of the physical processes that are crucial for 
understanding double compact object formation. 

In this paper we investigate the physics associated with 
the CE phase, focusing on the binding energy parameter, 
A, and the structure of the donor star. We introduce a 
new treatment of the NS/BH formation processes in core 
collapse supernovae, and assess the role of natal kicks, the 
maximum NS mass, and wind mass loss on the rates. To 
show the basic dependence of our results on metallicity, 
we consider two fiducial stellar populations: Z = Zq and 
Z = 0.1 Zq. Our main results are presented in Fig. [19] and 
[201 A broader metallicity study will be presented in paper 
II of this series, with inclusion of the cosmic star forma- 
tion history, and a calculation of the dependence of rates 
on redshift. Paper III will examine the detection progno- 
sis with advanced gravitational telescopes (e.g. Advanced 
LIGO/Virgo, Einstein Telescope), incorporating the lat- 
est gravitational waveforms and sensitivity curves for the 
various instruments. 

In order to make the results of our research eas- 
ily accessible to the community, we have provided 
an on-line Synthetic Universe database, available at 
www.syntheticuniverse.org. At present this contains 
detailed information on the population synthesis models in 
this paper. In addition, a number of physical parameters 
and their distributions that are not discussed here in detail 
(e.g., initial binary parameters, component masses, mass 
ratios, evolutionary tracks) are also provided, for those 
who might want to include them in their own research. We 
intend this database to act as an extension to our three 
papers, functioning as a repository of the full spectrum of 
our double compact object populations. The "Synthetic 
Universe" site will be updated in real time as new phys- 
ical models are calculated, and will respond to requests 
from the scientific community. Eventually, the Synthetic 
Universe will include various objects in addition to dou- 
ble compact object binaries, for example stellar popula- 
tions with white dwarfs (e.g., supernovae la progenitors; 
low frequency gravitational wave sources; cataclysmic vari- 
ables) and various types of X-ray binaries (e.g., persistent 
and transient sources; high- and low-mass X-ray binaries 
in Galactic and extra-galactic environments). 

In the next section we describe recent physical develop- 
ments addressing the CE events, the wind mass loss rates 
from massive stars, supernovae explosions and compact 
object formation, incorporated into the Star Track code. 
Section[3]contains a general description of the most impor- 
tant physical parameters used in our population synthesis, 
while in Section[3|we present the models used in this study, 
and our results. We summarize and conclude in Section[5] 

2. STARTRACK IMPROVEMENTS 

2.1. The StarTrack code 

We have used the population synthesis code StarTrack 
to calculate the numbers and properties o fDCOs. The full 
descri p tion o f the code can be found in iBelczvnski et al.l 
20021 120081). The code utilizes a set of stellar models 



Hurlev et al.l (|2000[ ); slightly modified from its original 



changes are described in Section 12.41 During core col- 
laps e, fallback and dir e ct BH formation is now accounted 
for (jFrver fc Kalogera . 2001) and the newly born objects 
receive natal kicks (jHobbs et al.l . |2005() . The formation of 
low mass NSs throu gh electron capture s u perno vae is also 
accounted for (e.g.. iPodsiadlowski et al.l (|2004[) ). Binary 
interactions are treated in detail, and the various physical 
processes have been calibrated using either resul ts of de- 
tailed evolutionary calculations (e.g., [Wellstein fc Langeij 
(I1999D for mass transfer sequences) , or specific sets of ob- 
servations fe.g.. lLevine et all (|2000[ ) for tidal interactions). 

In the following sections we describe the major updates 
to the code that are relevant for the formation of BHs 
and NSs. These crucial improvements are based on ad- 
vances in stellar evolution physics in the last decade. In- 
cluded are revised rates for wind mass loss for massive 
stars that now allow for the formation of BHs with masses 
observed in the Milky Way and other galaxies (see Section 
12. 2p . Another improvement relates to the CE coefficient, 
A, describing the binding energy of the envelope. It is 
now dependent upon the parameters (mass, radius, etc.) 
and the evolutionary stage of the donor. This is a ma- 
jor improvement when compared to most previous studies, 
where A was considered constant throughout the evolution 
of the donor (Section [^3)) . The third major update regards 
the core collapse/supernova explosions, and the resulting 
model for compact object formation. We introduce two 
convection-enhanced, neutrino driven supernova engines 
(based on Ray leigh- Taylor and Standing Accretion Shock 
instabilities) into our population synthesis code. This ad- 
dition allows us to account for the dearth of Galactic com- 
pact objects within the mass range 2-5 Mq (the so called 
"mass gap" fSection l2.4p ). 

2.2. Stellar winds 

The StarTrack code was recently updat ed with new wind 
mass lo ss rates (described in detail in IBelczvnski et al.l 
(l2010aD '). One of the major changes was the implemen- 
tation of mass loss for (D/B stars based on IVink et all 
(j200ll). followed by n ew w in ds for Wolf-Ray et star s 
(jHamann fc Koesterkel 119981; IVink fc de Koterl . I2005D . 
Fu rthermore, the Luminou s Blue Variable winds (based 
on iHumphrevs fc Davidson! (|1994[ )) were calibrated so as 
to allow for the formation of BHs with masses up to 15 Mq 
in a solar metallicity environment. This was done in order 
to account for th e observed masses of known BHs in the 
Milky Way (e.g. [O rosz et all (|201lD ). The most impor- 
tant consequence of these updates is the possibility of the 
formation of BHs with masses up to 30 Mq in sub-solar 
metallicity environments {Z = 0.3Zq), which corresponds 
to the observed mass of the most massive known stellar- 
mass BH in the ICIO X-1 bin ary (jPrestwich et al.l 120071; 
[Silverman fc FilipDenkol . l2008D . These revisions also result 
in the formation of BHs with masses up to 80 Mq in a low 
metallicity (Z = O.OIZq) environment. 

2.3. Common envelope updates 

We begin wit h an expression describing the energy bal- 
ance of the CE (IWebbinkl . [T98l : 



version) that allow for the evolution of stars at differ- 
ent metallicitics. The model for compact object forma- 
tion adopted in the code has been significantly revised; 
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Here, i?don.iob is the Roche lobe radius of the donor at the 
onset of Roche lobe overflow (RLOF), A is the binary sep- 
aration, Mdon.env is the mass of the ejected envelope, and 
M(jon and Mace are the masses of the donor and accretor, 
respectively. Indices i and / stand for initial and final pa- 
rameters, respectively. The parameter acE describes the 
efficiency of the transfer of orbital energy into unbinding 
the envelope. If there is sufficient orbital energy to un- 
bind the envelope the binary may survive the CE phase, 
otherwise the component merger and the formation of a 
peculiar single object is assumed (an example: if the ac- 
cretor is a compact object then in the outcome of the CE 
a star harbouring a NS /BH may form) . The value of ac e 
has been set to 1 throughout this study. The range of val- 
ues for the product acEA may therefore be covered by A 
alone. The parameter A describes the binding energy of 
the envelope, and is defined as: 

J^bind = , i^j 

where R is the radius of the donor, which becomes equal to 
the Roche lobe radius of the donor at the onset of the CE. 
The most significant changes of binding energy are usu- 
ally associated with radial expansion or contraction of a 
star. Additionally, the mass loss/gain may alter the bind- 
ing energy. Generally, as the envelope of the star expands 
it becomes less dense and so its binding energy decreases, 
while a contraction of the envelope leads to an increase 
in the binding energy. The parameter A is a function of 
binding energy, radius, and mass (Eq. [2]). For example, A 
will decrease if the radius expands faster than the bind- 
ing energy decreases , which is the case for stars on the 
Hertzsprung gap (HG) stars (Fig. [THll)- If the radial ex- 
pansion of a star remains moderate, A may increase as it 
happens on the Red Giant Branch (RGB) and Asymptotic 
Giant Branch (AGB) (see Fig. [Ij. For the most massive 
stars (initial mass above ^ 20 M0) the mass loss becomes 
a significant factor, and the changes in density profile, ra- 
dius, and binding energy result in an almost constant A 
during later evolutionary stages (see Fig.[2H4]). 

Most importantly, Eq. [2] implies that given the mass, 
radius, and evolutionary state (i.e. the envelope size) of a 
donor, a change of A will change the binding energy of its 
envelope. For example, the binding energy of the envelope 
of a star with A = 0.1 will be larger by an order of magni- 
tude as compared to an identical star with A = 1. During 
the CE phase, orbital energy is used to eject the envelope. 
Therefore, a binary will have much more difficulty eject- 
ing the tightly bound envelope of a donor with A = 0.1, as 
compared with one with A = 1. 

The consequences for the formation of double compact 
objects are twofold. If A has a low value, and the stars have 
insufficient separation (orbital energy), they will merge 
during the CE phase. This terminates further binary evo- 
lution and prevents the formation of a DCO. If, on the 
other hand, A has a high value, it is most likely that a given 
system will eject the CE with a small resulting decrease 
in orbital separation. Binaries (later becoming DCOs) re- 
taining a wide enough separation after the CE event may 
not be able to merge within a a Hubble time. Such DCOs 
become uninteresting from the gravitational wave detec- 
tion perspective. Also, a supernova explosion occurring in 



a binary that maintains a wide orbital separation tends to 
disrupt the system. 

2.3.1. Hertzsprung gap donors 

The outcome of the CE phase depends strongly on 
th e evolutionary stage of the donor, as first discussed 
in iBelczvnski et al] (|2007[ ). Even before considering the 
energy balance, it is critical to incorporate an understand- 
ing of the core-envelope structure of the star. For example. 
Main Sequence (MS) stars do not have clear core-envelope 
division, as the helium core is still in the process of be- 
ing developed. Stars on the HG similarly lack a clear 
entropy jump associated with the core-envelope structure 
(jivanova fc Taanil . |2004|) , although when in the evolution 
such a division appears (late HG or post HG?) remains 
unclear. In the case of a CE initiated by a MS or HG 
donor, therefore, the orbital energy is transferred into the 
entire star rather than just the envelope. This makes the 
ejection of the envelope difficult, and for MS stars we as- 
sume that this will always result in a merger. However, 
in the case of HG we extend the analysis by consider- 
ing two possibilities. One is to ignore the core-envelope 
boundary issue, and proceed with the calculation of the 
energy balance (submodel A). The second is to adopt a 
pessimistic approach (submodel B) in which each CE with 
a HG donor leads (independent of energy balance) to a 
merger. Both submodels are calculated for each evolu- 
tionary model in this study. Massive stars, beyond HG, 
that enter the core helium burning already possess well- 
defined core-envelope structure and require only the en- 
ergy balance calculation. Exceptions are the Helium MS 
and Helium HG stars, which have an analogous structure 
to their hydrogen counterparts, and are treated as such. 

2.3.2. New calculation of the A parameter 

A physica l estimate of A was recently presented b y 
two groups: IXu fc Lil ()2010D and iLoveridge et al.l (|2011[ ). 
We adopt the expressions from the former group in the 
StarTrack code, although both calculations yield similar 
results for NS/BH progenitors. The new A values cover 
CE events, and depend on the evolutionary stage of the 
donor, its mass at Zero Age Main Sequence (ZAMS) and 
the mass of its envelope, and its radius. In addition, all of 
these quantities can depend on a wide range of metallicities 
{Z = 10^^-0.03). This represents a significant improve- 
ment in the physics of the formation of DCOs, as most 
previous studies treated A as constant (usually A = 1) 
throughout the evolution of a star. In the new approach 
the authors calculate two values of A from detailed stellar 
models: in the first they treat the binding energy of the 
donor's envelope as consisting only of its gravitational en- 
ergy, while in the second the binding energy is decreased 
by the full internal energy of the envelope. S ince neither 
of these extremes is plausible (|Xu &: Lil . |2010() , we choose 
our A to be the average of the two. In the case of stars 
with masses leading to the formation of BHs we find that 
A <J 1 (see Fig. I3l4p . We find higher values for typical NS 
progenitors, and especially in the late evolutionary stages 
we have A ^ 1 (see Fig. [T]). Since the group that devel- 
oped this procedure is affiliated with Nanjing University, 
we label this approach as the Nanjing A. We have addi- 
tionally requested that the authors compute a few extra 
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models for massiv e stars, as thes e arc particularly rele- 
vant to our study. IXu fc Lil (|201Cl[ ) presents models up to 
ams — 20 Mq, and we have now obtained models up to 
-^zams = 100 Mq. We then extrapolate these to our en- 
tire mass range (up to Mzams = 150 Mq). For stars with 
Afzams > 100 Mq we usc the fitting function calculated 
for Mzams = 100 Mq, but input the higher mass when 
necessary. 

2.4. Compact object formation 

To evolve a star from ZAMS toward its eventual super- 
nova (SN) wit h the StarTrack code, w e use the (slightly 



modified; see Beiczvnski et al.l (|2002| )') procedure pre- 
sented in iHurlev et al.l (|2000[ ). with updated wind mass 
loss rates as described in Section 12.21 The mass of the 



compact object is calculated from the properties of the 
pre-supernova star, and the type (whether it is a NS or a 

BH) is set solely by its mass^ 

We use a recent study bv lFrver et al.l (j201ll ) to describe 
the SN explosion and the resulting compact object forma- 
tion. Our models allow for a successful explosion without 
the need for the artificial injection of energy into the ex- 
ploding star. This is a major update of the input physics 
used in the population synthesis of massive binaries. Pre- 
vious studies at best used significantly outdated supernova 
models, if they used any at all. We have introduced two al- 
ternative supernova models into our code: "Delayed" and 
"Rapid" . Both are core-collapse scenarios, and they share 
the same convection-enhanced neutrino-driven explosion 
mechanism. The main factor that differentiates the two 
models is the type of instability which causes the macro- 
scopic flows of matter (similar to convection) that eventu- 
ally lead to the ejection of the infalling matter. The De- 
layed model is sourced from the standing accretion shock 
instability (SASI), and can produce an explosion as late 
as 1 s after bounce, while the Rapid model starts from the 
Ray leigh- Taylor instability and occurs within the first 0.1- 
0.2 s. In the Rapid scenario we either end up with a very 
strong (high velocity kick) supernova in the case of a low 
mass star (Mzams ^ 25 Mq), and produce a NS, or the 
supernova fails, and there is direct collapse to a massive 
BH. In the Delayed case the entire spectrum of explosion 
energies is allowed, and this results in a wide range of 
compact object masses, from NSs to light BHs to massive 
BHs. The formulae describing the two supernova engines 
adopted in StarTrack arc provided in iFrver et al.l (j201lh 
(Eqs. 15-20). 

We also allow for the form ation of NSs through Electron 
Capture Supernovae (ECS. iMivaii et al.l (|l980f )V These 
are weak supernovae (no natal kick assumed) occurring 
for the lowest mass stars (Mzams ~ 7 Mq), and they end 
up forming NSs. 

3. MODELING 

The physics underlying the formation of double compact 
objects is uncertain, notably due to modelling challenges 
of the supernova and CE phases. Core collapse, usually 
(but not always) followed by a supernova explosion, forms 
a compact object. The mass and the type of compact 
object is determined by the details of the event. Addi- 
tionally, supernova asymmetry via mass loss and/or a na- 
tal kick may disrupt a binary, depleting the population of 



double compact objects. At the same time, for some bina- 
ries a natal kick of the right amplitude and direction may 
produce a prematurely coalescing DCO. The CE phase is 
present in all formation scenarios, for all types of close 
double compact objects. It is the primary mechanism for 
bringing initially widely-separated binaries into close or- 
bits, thereby allowing them to coalesce within a Hubble 
time. 

To explore the uncertainties associated with our models 
we calculate a suite of population synthesis results, in- 
vestigating a range of factors which have the largest im- 
pact on the rates and physical properties of DCOs. As 
a reference we use the standard model introduced in Sec- 
ti onl4.1[ that resembles i nput physics described in detail 
in lBelczvnski et all (j2008| ). with some important additions 
and modifications. The subsequent models are variations 
on this standard model, each exploring a single parameter 
connected to cither the core collapse/supernova explosion 
or CE phases. Despite the fact that we now have phys- 
ically motivated values for A, we also explore models in 
which A is constant, with the specific value set over a wide 
range. The range is chosen to be such as to encompass all 
plausible values allowed by detai l ed ph ysical calculation 
of the Nanjing results ijXu fc Lil . |2010[ ). Variations 1-4 
employ A =0.01, 0.1, 1, and 10, respectively (see TablejT]). 

To delineate between a NS and a BH we need to adopt 
a value for the maximum NS mass. Theoretical studies of 
the equation of state a Uow v alues in the range 1.5-3.0 Mq 
(jLattimer fc Prakashl |2010[ ). Observations yield a nar- 
rower range, with the most massive NSs reaching 2 Mq 
(|Demorest et al.l . |2010[) . Statistical analyses of the mea- 
sured BH masses indicate that the BH mas s distribution 
is unlikely to extend below about 4.5 Mq ijBailvn et al.l 
Il998t lOzel et all . 120101: iFarr et all . l201lt) . which leaves a 
significant range with no compact objects. This might 
argue for a higher upper limit on the NS mass (or even 
potentially a lower limit on the BH mass?). Utilizing the 
theoretical and observational estimates, we vary the maxi- 
mum NS mass from 3.0 Mq (Variation 5), through 2.5 Mq 
(Standard); and as low as 2.0 Mq (Variation 6). 

In the standard model for natal kicks in core collapse su- 
pernovae we employ a Maxwellian kick distribution with 
a = 265 km s~^ , based on observe d velocities of single 
Galactic pulsars (|Hobbs et al.l . I2005D . The magnitude of 
the kicks is modified by the amount of fallback during a 
SN as follows: 



Vfe — Vmaxi^ — ffb), 



(3) 



where Vk is the final magnitude of the natal kick, Vmax 
is the velocity drawn from a Maxwellian kick distribution, 
and ffb is the fallback factor. The values of ffb range be- 
tween 0-1, with indicating no fallback/full kick and 1 
representing total fallback/no kick (in this case a "silent 
supernova" ) . This factor is calcula ted for both t he De layed 
and Rapid SN engines according to lFrver et al.l (|2011l ). For 
the electron capt ure supernovae the e xplosions are found 
to be symmetric (jDessart et al.l . |2006() , and we assume no 
natal kick. However, the orbit is still modified due to the 
mass loss in the explosion. BHs that form through partial 
fall back are assumed to receive natal kicks that arc drawn 
from the same distribution as for the NSs, but the value is 
decreased in proportion to the amount of fall back. In par- 
ticular, the most massive BHs form via full fall back, which 
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leads to no natal kick and "silent" B H formation. This 
may be supported both the oretically (iFrver &: Kalogerai 
l200lh and from observations ([Mirabel fc Rodriguesl . l2003| ). 
There seems to be some evidence that NSs that 
fo und in binaries receive smaller nat al kick s 
Belczynski et alj (|2010c[ ): IWong et all (|2010f ): 



are 
(e^ 



iBodaghee et all (|201lt)l. of the order of 100 km s~^ Keep- 



ing the same Maxwellian distribution as in the standard 
model , we modify the magnitude of natal kicks of NSs 
and BHs to cr = 132.5 km s~^ (Variation 7). 

Additionally, we calculate two models in which the BH 
kicks take on two extreme values. In one, newly formed 
BHs receive full natal kicks as, observed for single pulsars 
(cr = 265 km s^^) independent of the amount of fallback 
(Variation 8). In the other, no natal kicks are applied to 
BHs, independent of their mass or the amount of fallback 
(Variation 9). 

As our standard model we select the Rapid supernova 
engine, since this supernova mechanism, combined with 
binary evolution, succes sfully reproduces the mass gap 
(|Belczvnski et al.l l2011bl ) observed in Galactic X-ray bi- 
naries (|Bailvn et all . Il998t lOzel et all . l2010t) . On the 
other hand, the Delayed engine (Variation 10), as well 
as the routine previous ly implemented within Star Track 
(jBelczvnski et al.l . l2002[ ). generates a continuous spectrum 
of compact object masses. 

Observational and theoretical developments have led 
to the discovery of a number of puzzling phe- 
nomena in the winds from massive stars. The 
two most prominent are the "wea k wind problem" 
(e.g., IChletjowski fc Garmanvl (|1991[) : iKudritzki et al.l 



19911): iHerrero et al.l (120021)') a nd "wind clu mping" 



C^g.-lOsterbrock fc Flatheil (ll959D:lMarkova et al.l (12004... 
iRepolust et al.l (|2004[) : iLepine fc Moffat! (|2008D ). The for- 
mer is related to the fact that wind mass loss rates from 
late O and early B type stars might be ^ 1-2 orders of 
magnitude lower than theoretically predicted. The latter 
suggests that winds might be forming dense clumps rather 
than being distributed uniformly, which may lead to an 
overestimate of mass loss rates by a factor of <; 2. Based 
on this we investigate the possibility that our wind mass 
loss rates are too high, and so we calculate a model (Vari- 
ation 11) in which the rates are reduced by a factor of 
2. This is done for all stars at all points in their nuclear 
evolution. 

Table [T] lists all of our models, each with the relevant pa- 
rameter to be varied as compared to our standard model 
(listed in the first row). Each model is calculated for two 
metallicity values: solar {Z = = 0.02) and 10% so- 
lar {Z = 0.1 Zq 0.002). Each model is further divided 
into submodel A (CE energy balance with HG donor) and 
submodel B (CE merger with HG donor) as discussed in 
Section [2.3.11 We end up with 45 distinctive models (11 
variations for 2 metallicities and 2 CE submodels, on top 
of our standard model). 

For each model we evolve 2 x 10^ binaries (with one ex- 
ception, see. Sec. 14. 5p . assuming that each component is 
created at the same time. Each binary system is initialized 
by four parameters which are assumed to be independent. 
These are: primary mass Mi (initially more massive com- 
ponent), mass ratio q = M2/M1, where M2 is the mass of 
the secondary component (initially less massive), the semi- 



major axis a of the orbit, and the eccentricity e. The mass 
of the primary component is rand omly chosen from th e 
initi al mass function adopted from iKroupa et al.l (|1993D . 
and lKroupa fc WeidneH (|2003D . 
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0.08 M© < Ml < 0.5 Mq 
0.5 Mq < Ml < 1.0 Mq 
1.0 M© < Ml < 150 Mq, 



(4) 

where a = 2.7 is our standard choice for field popula- 
tions. Stars are generated from within an initial mass 
range Mmm^Mmax, with the limits based on the targeted 
stellar population. For example, NS studies require evo- 
lution of single stars within the range 8-25 Mq, while 
for BHs the lower limit is 25 Mq. Binary evolution may 
broaden these ranges due to mass transfer episodes, and we 
therefore set the minimum mass of the primary to 5 Mq. 
We assume a flat mass ratio distribution, ^{q) = 1, over 
the ra nge q = 0-1, in agreeme nt with recent observa- 
tions (jKobulnickv fc Frveij . l2007f ). Given a value of the 
primary mass and the mass ratio, we obtain the mass of 
the secondary from M2 = qMi. However, for the same rea- 
sons as for the primary, we don't consider binaries where 
the mass of the secondary is below 3 Mq. The distribu- 
tion of initial bina ry separations is assumed to be flat in 
log(a) (lAbtl . ll983D . and so cx -, with a ranging from val- 
ues such that at ZAMS the primary fills no more than 50% 
of its Roche lobe to 10^ Rq. For the initial ecce ntricity 
we ad o pt a thermal equilibrium dist ribution fe.g.. [Heggid 
dHli); iDuquennov fc Mavod (|l99lf )) S(e) = 2e, with e 
ranging from to 1. 

4. RESULTS 

The double compact object merger rates for a synthetic 
galaxy resembling the Milky Way, but for two differing 
values of metallicity, are presented in Table [2] (Z = Zq) 
and Table|3|(Z = 0.1 Zq). The corresponding plots of the 
merger rates are shown in Figures [T9| and |20l The physi- 
cal properties of double compact objects for the standard 
model are presented in Figure|7|(chirp mass), Figure|5|(de- 
lay time), and Figure |9| (mass ratio). The distribution of 
chirp masses (Figures [TOHTSi (with corresponding Tables 
IIHH]) and delay times (Figures [T4lfT7|) for all the models 
are presented. The delay time, tdei, is the sum of the time 
needed to form two compact objects from a ZAMS bi- 
nary and the time for the two compact objects to coalesce 
due to the emission of gravitational radiation. For double 
compact objects the former evolutionary time interval (~ 
Myr) is usually much shorter than the latter merger time 
(~ Gyr), and the delay time is rather similar to the merger 
time. Additional models and DCO population properties 
are available on-line at www.syntheticuniverse.com. 

For each model we calculate the Galactic merger rates. 
These are defined as the number of coalescences of DCOs 
per unit time occurring in a synthetic galaxy similar to the 
Milky Way (with age of 10 Gyr and a constant star forma- 
tion rate (SFR) of 3.5 Mq yr~-^). In practice this is done 
by checking if the delay time of a DCO (with a random 
starting point between 0-10 Gyr) falls near the current 
galaxy age (10 Gyr). However, the amount of mass within 
the simulated binaries corresponds only to a part of the 
star forming mass of the galaxy. In order to extrapolate 
the simulated mass to that of the entire galaxy we em- 
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ploy the following procedure. First, the amount of mass 
contained within the simulated binary stellar population 
is estimated (the mass of the primary follows from Eq. |4l 
and the mass of the secondary follows from our assumed 
mass ratio distribution). Additionally, we assume a bi- 
nary fraction of 50%, so that for each binary system there 
is one additional individual star. The mass of each of the 
individual stars is taken to be in the same range as for pri- 
mary components in binaries. The acquired mass, Afacq, 
is then divided by the age of the synthetic galaxy tgai (10 
Gyr) in order to get a constant star formation rate corre- 
sponding to the simulated stellar mass (0.073 Mq yr~-^). 
To match this SFR to the one of the synthetic galaxy 
SFRgai (3.5 M0 yr~^) one needs a multiplication factor 
of 48 (/sfr)- The corresponding equation is: 

/sFR = SFRgai (^^) ' (5) 

Therefore, to extrapolate our results to the entire mass in 
the synthetic Milky Way, we use each synthetic DCO bi- 
nary 48 times. Each time the given binary is assigned a 
new starting time (from a uniform distribution), and if its 
coalescence time falls within 9-10 Gyr it is included in our 
results. 

The typical range of the number of DCOs, N, gener- 
ated in each simulation is ~ 1000-10000 for submodel A 
and ^ 10-100 for submodel B. Therefore the relative sta- 
tistical error (v^ /N) is at most '--^ 10%-30% where the 
ranges are given by the values in submodels A-B). The 
errors arising from uncertainties in various aspects of the 
single and binary star evolution can change the Galactic 
merger rates of DCOs by ^ 1-2 orders of magnitude (as 
shown in the next sections), making the statistical errors 
irrelevant. 

In addition to the population of DCOs with delay times 
below 10 Gyr (the merging population), in each model we 
acquire another population occupying the domain above 
this time limit (the non-merging population). This pop- 
ulation contains each type of DCO, and is available at 
www . syntheticuni verse . org. 

4.1. Standard Model 

At solar metallicity, the merger rates for a synthetic 
galaxy similar to the Milky Way are dominated by NS- 
NS systems (23.5-7.6 Myr'^; submodel A-B), with a 
smaller but still significant contribution from BH-BH sys- 
tems (8.2-1.9 Myr-i), and with a minor contribution from 
BH-NS systems to the overall DCO merger rate (1.6 — 0.2 
Myr"-*^; see Tabic [J). Qualit atively these findings ar e con- 
sistent with previous results ([Belczvnski et al.l |2002[) . The 
quantitative results, however, are quite different, due to 
the many improvement s in th e mod els over the interven- 
ing decade. iBelczvnski et al.l ()2002[ ) found the following 
mean rates: NS-NS 53 Myr"!, BH-BH 26 Myr"! and 
BH-NS 8.1 Myr"-'^. At sub-solar metallicity, the systems 
with BHs increase their relative contribution to the overall 
rates, and the merger rate is dominated by BH-BH systems 
(73.3-13.6 Myr~^), with smaller contribution from NS-NS 
(8.1-2.5 Myr-i) and BH-NS systems (3.4-2.3 Myr'^; sec 

^This factor is the ratio of the number of stars between 8-20 Mq 
progenitors) as calculated from the IMF (see Section [Sj. 



Table [3]). These results arc qualitatively similar to our re- 
cent work on the depende nce of merger rates on metallicity 
(|Belczvnski et al.l 120101 /) . Agai n, quantitativel y there are 
significant differences, feelczvns ki et al.l (j2010b[ ) found the 
following rates: 84-6.1 Myr^^ (BH-BH), 41-3.3 Myr'^ 
(NS-NS) and 12-7.0 Myr"! (BH-NS) for Z = 0.1 Zq. 
These changes reflect the fact that since the previous study 
we have introduced physical A values and observationally 
constrained SN models, which yields a new compact object 
mass spectrum. 

A general division of these rates by DCO type may be 
understood in the following way. The initial mass function 
(IMF) falls steeply with mass, and delivers more NS than 
BH progenitors (by a factor Q of ^ 4) . The supernova 
explosion is the major process that drastically affects the 
number of massive binaries , as the explosions tend to dis- 
rupt binaries. This is especially true in the case of NS pro- 
genitors, as these receive large natal kicks and as many as 

90-95% of potential binaries may end up disrupt e d afte r 
the first supernova expl osion (e.g., iLorimer et al.l ()2004[ ): 
IBelczvnski et al.l ()2010cl )). Supernovae do not affect bina- 
ries with BHs as much because it is believed that most of 
the massive stars producing BHs do not experience large 
kicks at core collapse. What follows is that the rates for 
NS-NS and BH-BH mergers are not as separated as would 
be simply deduced from the IMF. The rates are also very 
sensitive to the details of the CE, and these are quite dif- 
ferent for NS and BH progenitors (as discussed below). 
Additionally, the three types of double compact objects 
evolve along separate evolutionary channels (see Table |4]) . 
This qualitative picture can explain the calculated ratio 
of NS-NS to BH-BH merger rates (~ 3-4). The BH-NS 
merger rates are the smallest, as the majority of potential 
progenitor binaries have initially large mass ratios, and 
the first interaction of the two components leads to a CE 
phase and the inevitable merger. This happens because 
the massive envelope of the BH progenitor cannot be suc- 
cessfully dispersed by the much less massive NS progenitor. 
Note that none of the significant BH-NS formation chan- 
nels starts with the CE phase, but instead the progenitors 
of these systems originate from a narrow mass range (mass 
ratio of the two components larger than ^ i-i allowing 
for the first interaction to be a stable RLOF (see Tabled]). 

At sub-solar metallicity other factors come into play and 
make BH-BH systems dominant in the overall merger rate. 
First, the smaller wind mass loss makes prc-SN progeni- 
tors of NSs slightly more massive (by ^ 10%). Heavier NS 
progenitors tend to explode as core collapse supernovae 
(full kicks) rather than ECSe (no kicks, sec Section [2^ . 
This means that for sub-solar metallicity more NSs are 
formed with disruptive natal kicks than for a solar envi- 
ronment. Additionally, smaller wind mass loss decreases 
the expansion of the separation between the components. 
This causes the progenitors of NS-NS systems, at sub-solar 
metallicities, to engage in a second CE (see Table [5]) just 
after the first one. This increases the probability of a 
merger during evolution when compared with Zq, hence 
the drop of merger rates of NS-NS systems from Zq to 
0.1 Zq. The merger rates for BH-BH systems increase for 
sub-solar metallicity, as low wind mass loss rates allow the 

(NS progenitors) to the number of stars between 20-150 Mq (BH 
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progenitors to remain more massive during evolution. This 
in turn makes the pre-SN stars more massive and allows 
for larger amount of fallback. Increased fallback reduces 
the magnitude of the natal kicks (sec Eq. |3]) to almost none 
or none at all (direct collapse into a BH) and makes the 
SN significantly l ess disr uptive, as explained in detail by 
iBelczvnski et all (|2010b[) . 

The formation of BH-NS binaries is determined by the 
properties of the progenitors of both compact object types. 
Therefore, the behaviour of these systems may be consid- 
ered as a combination of effects noted in the formation of 
BH-BH and NS-NS systems. 

In Figures [SJ [51 [71 HI and [HI we present the distributions 
of DCO progenitor ZAMS masses, DCO masses, chirp 
masses, delay times and DCO mass ratios for the stan- 
dard model. For purposes of illustration, the distribu- 
tion of progenitor and remnant masses are given for sub- 
model A, Z = Z0 only. Additional plots are available at 
www. syntheticuni verse . org. 

In the distribution of progenitor ZAMS masses (Fig. [S|) 
one can clearly see that binary evolution blurs the limits 
for ZAMS mass of the star for the formation of NS/BH. 
For NS-NS progenitors (top panel), the masses of the 
primary components (up to ~ 30 Mq) exceed the typi- 
cal upper limit for the formation of a NS for single stars 
('^ 20 M0 , also lower limit for the formation of a BH) . Un- 
der favourable circumstances, binary evol ution may push 
this limit even further, up to 100 M ^ (e.g.. lWellstein etaTI 
(|200lh : IBelczvnski fc TaamI (|2008[) ). Note that for BH- 
NS systems (middle panel), a high progenitor mass does 
not necessarily imply that it will form a BH. During their 
evolution, progenitors of BH-NS systems may undergo a 
mass ratio reversal (due to mass transfer events) so the 
primary component (initially more massive) may become 
a NS (as seen on the middle panel of Fig. [6]). For BH- 
BH progenitors the lowest mass of the primary compo- 
nent is ^ 45 Mq , and the lowest mass of the secondary is 
~ 25 M0. Also, progenitor stars of BH-BH systems have a 
wide mass distribution (ranging up to 150 Mq). However, 
the masses peak at ~ 55 Mq for the primary component 
and ~ 40 Mq for the secondary. 

The distribution of masses of DCO remnants (Fig. [5]) 
clearly shows the gap between the upper mass achiev- 
able by NSs (2 Mq, despite allowing for the formation of 
NSs with mass up to 2.5 Mq) and the lower mass of BHs 
(~ 5 Mq). This is due to the implementation of the Rapid 
supernova engine, which is the first to successfully repro- 
duce this "mass gap" (for details see Section 12. 4p . Ad- 
ditionally, for BH-NS binaries some remnants formed out 
of the initially less massive star in the binary fall within 
the BH mass regime while some primary components (ini- 
tially more massive) may end up as NSs. This is due to 
the aforementioned mass ratio reversal. Despite a wide 
range of BH-BH progenitor masses, the component masses 
of remnant systems are mostly clustering around 5-9 Mq . 
Such a drastic reduction in mass range comes from the 
significant wind mass loss for massive BH progenitors and 
mass ejection in CE events. Both factors usually reduce 
the masses of the whole progenitor stars to the masses of 
their cores, for which the mass range is narrow. 

To calculate the chirp mass Afchirp of a DCO we use the 



following formula: 

Mehhp = (MiM2)t (Ml + Ahyi, (6) 

where Mj are the masses of the components. The most 
notable aspect of the chirp mass distributions (Fig. [7]) is 
the maximum value for BH-BH systems (~ 30 Mq), for 
submodel A, 0.1 Zq. As metallicity decreases from solar 
to sub-solar, so do the wind mass loss rates. This allows 
for the formation of more massive BHs. For submodel B, 
Z — 0.1 Zq the maximum value is lower (18 Mq) than 
for submodel A for the same metallicity. In submodel A 
the most massive progenitors of BHs experience significant 
expansion during evolution leading to the CE events with 
a HG donor, rather than a donor beyond the HG which 
is more likely for lower mass progenitors. Submodel B 
does not allow for the survival of HG donors during a CE 
(see Section [2.3.1^ and the most massive BHs disappear. 
The average chirp masses for Z ~ Zq are the same for 
both submodels: 1.05 for NS-NS, 3.2 Mq for BH-NS, and 
6.7 Mq for BH-BH systems. For 0.1 Zq the average chirp 
masses are: 1.1-1.1 Mq (submodel A-B) for NS-NS, 3.2- 
3.1 Mq for BH-NS, and 13.2-9.7 Mq for BH-BH systems. 

The delay time distributions (Fig. [5]) are proportional to 
f^^^. The average delay time for systems merging within 
10 Gyr at Zq is: 1.1-1.5 Gyr (submodel A-B) for NS-NS, 
1.7-1.7 Gyr for BH-NS, and 1.0-2.5 Gyr for BH-BH sys- 
tems. For sub-solar metallicity the average is: 1.0-2.3 Gyr 
for NS-NS, 1.5-1.5 Gyr for BH-NS, and 1.0-1.4 Gyr for 
BH-BH systems. 

The delay times for submodel B are in general higher 
than for submodel A. This is a direct result of our as- 
sumption on the outcome of the CE phase with HG donor. 
Binary populations are identical in submodel A and B un- 
til a CE event. If the separation is relatively small, the CE 
may be initiated early on in the evolution, specifically by a 
HG donor. The survival of such an event may lead to the 
formation of a close binary in submodel A. In submodel B 
CEs with this type of donor are always considered a merger 
(see Section [533J ■ Binaries with separation large enough 
to prevent the rapidly expanding HG star from overfill- 
ing its Roche lobe will initiate the CE in later stages of 
evolution. This scenario meets the criteria of A and B sub- 
models and allows the binaries to be accounted in both (in 
terms of forming merging DCOs). Submodel B, therefore, 
favours binaries with larger separations, which translate 
into larger merger times. 

Additionally, the delay times decrease with metallicity. 
This comes from the fact that for sub-solar environments 
stars lose less mass in winds and therefore, form more mas- 
sive remnants. The more massive the components of a 
DCO, the less time it takes for a system to merge. A 
secondary effect of the reduced wind mass loss rates for 
Z < Zq on the delay times is the smaller expansion of 
the separation during evolution as a smaller amount of 
mass is removed from the system. 

The distribution of mass ratios q is similar for both 
submodels and both metallicities. NS-NS systems group 
around q = 0.9 as the NS masses found in our models are 
on average similar in each formed DCO system {'^ 1.1- 
1.4 Mq; see Fig. [5]). BH-NS systems group around small 
mass ratio values {q — 0.2), which means a large difference 
in mass between both remnant compact objects. This fol- 
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lows from the fact that the typical NS mass is 1.3 Mq 
and the typical BH mass is found at 5-9 Mq (see Fig. [5]). 
BH-BH systems have the widest range of q, and typically 
our simulations show mass ratios in the range 0.4-1.0 with 
increasing probability toward q = I. This simply reflects 
the fact that BH progenitors come from a wide range of 
masses, and that binaries with similar masses more readily 
survive binary interactions (mass transfers and supernovae 
explosions). 

4.2. Variation 1 

This is the first of four variations addressing the CE 
binding energy parameter, A. In this model we fix A = 
0.01. This It is found that in this evolutionary variation 
all merger rates significantly decrease as compared to the 
standard model (see Table [D and [3]) . For solar metallic- 
ity the Galactic merger rates are dominated by BH-BH 
systems (1.1 Myr~^ for both submodels), followed by NS- 
NS (0.4 Myr"^ for both submodels), and BH-NS systems 
(0.002 Myr-i , also for submodels A and B). For 0.1 Z© 
the rates for BH-BH systems are 12.5-8.1 Myr~^ (sub- 
model A-B), 0.06 Myr-i for NS-NS systems (for both 
submodels), and 0.03 Myr"-'^ for BH-NS systems (also for 
both submodels). The behaviour of the population of dou- 
ble compact objects in this model can be understood as 
follows. We have chosen A to have a very low value, which 
translates into a very high binding energy of the CE (see 
Eq. [21). This binding is so strong that most of the bina- 
ries experience a merger during the first CE phase in their 
evolutionary history. This terminates the evolution of the 
binary and prevents the formation of a DCO. The least 
affected are BH-BH systems, since the fixed A value used 
in this variation is comparable to (although lower than) 
the Nanjing A (see Fig. I3l4p of the typical BH progenitor. 

In earlier stages, such as on the Hertzsprung gap, the 
values are higher [Nanjing A « 0.1-0.2; see Fig. [3]). This 
means that in this fixed A scenario the value of the binding 
energy of the CE of a HQ donor is an order of magnitude 
larger than for the standard (more physical) model. There- 
fore, many CE events with a HG donor end in a merger, 
as the reservoir of the orbital energy is insufficient to eject 
the envelope. The surviving binaries are the ones in which 
the CE was initiated by a CHeB donor and/or that have 
very massive accretors and therefore a large orbital energy 
reservoir. Submodel B does not allow for CE events with 
a HG donor. In effect, the very small A (very large bind- 
ing energy) adopted in submodel A also prevents binaries 
from surviving a CE event with a HG donor, just as in 
submodel B. As a result, the binary populations formed in 
the two submodels are very similar. 

For Z = Tjq the average chirp mass for NS-NS systems 
is 1.1 Mq, 3.2 Mq for BH-NS, and 6.5 Mg for BH-BH 
systems (for both submodels). In both submodels the dis- 
tributions of chirp masses in the standard model and VI 
are similar. The slight differences come from the fact that 
the population of merging DCOs in VI is composed of bi- 
naries that would form the non-merging population in the 
standard model. For Z = O.l Zq the corresponding values 
are 1.1 Mg for NS-NS and 3.6 M© for BH-NS systems (for 
both submodels). In the case of BH-BH systems the av- 
erage chirp mass is 20.0 Mg for submodel A and 16.1 Mg 
for submodel B. For Z = 0.1 Zg, submodel A, the dis- 
tribution of BH-BH chirp masses is flatter in VI than for 



the standard model. Specifically, there are fewer low chirp 
mass systems in VI than for the standard model. There 
are two mechanisms determining this outcome. The first 
is the increased number of mergers during the CE in this 
variation, which reduces the overall number of DCOs. The 
second is the fact that less massive BH-BH systems have 
relatively light progenitors. These binaries have smaller 
chances of ejecting such a strongly bound CE due to a 
smaller orbital energy reservoir, and merge in the process. 

Since the populations in submodels A and B are simi- 
lar, we find similar distributions of chirp masses and delay 
times. For Z = Zg the average chirp mass for NS-NS 
systems is 1.09 Mg, 3.2 Mg for BH-NS, and 6.5 Mg for 
BH-BH systems. For Z = 0.1 Zg the values are 1.11 Mg 
for NS-NS, 3.6 M© for BH-NS systems for both submod- 
els, while for the BH-BH systems the average chirp mass 
is 20.0 Mg for submodel A and 16.1 Mg for submodel B. 
The components of these BH-BH systems belong to the 
most massive compact objects, and this results in the high 
cutoff point of the chirp mass distribution in submodel B 
for 0.1 Zg (~ 30 Mg), when compared to other variations. 
Due to the very low A in VI, only very wide progenitor bi- 
naries can survive a CE phase. These binaries would form 
non-merging (tdci > 10 Gyr) BH-BH systems in the stan- 
dard model. Since the progenitor binaries are very wide, 
the CE is typically initiated when the donor expands sig- 
nificantly. This happens at late evolutionary stages and/or 
for very massive donors; therefore, the donor in most cases 
has already evolved past the HG and is in the CHeB phase. 
As most donors are CHeB stars, both submodels produce 
very similar populations. Furthermore, since many pro- 
genitors are very massive stars, they produce very massive 
BH-BH systems in both submodels. 

The average delay times for Zg are ^ 52 Myr for NS-NS, 
- 1.7 Gyr for BH-NS, and ~ 1.5 Gyr for BH-BH systems. 
For 0.1 Zg these arc - 82 Myr for NS-NS, - 0.8 Gyr for 
BH-NS for both submodels. For BH-BH systems we find 
average delay ~ 480 Myr for submodel A and ~ 680 Myr 
for submodel B. The short delay times for NS-NS systems 
is a direct consequence of the very low value of A. The high 
binding energy of the envelope following from this causes 
significant orbital energy dissipation. This prevents the 
formation of NS-NS systems with delay times over 300 
Myr (sec Fig. [HHITI). 

Due to the severe reduction in the total number of sys- 
tems in this model, the least populated group, BH-NS sys- 
tems, is subject to larger statistics errors. For example, for 
solar metallicity we find only 1 merging BH-NS system in 
the VI model. However, as we will see below, this model 
can be excluded based on observations of known NS-NS 
Galactic systems, and therefore we do not follow up with 
additional computations. 

4.3. Variation 2 

In this model we fix A = 0.1. The rates for solar metal- 
licity, for the DCOs are: 11.8-1.1 Myr^^ for NS-NS, 2.4- 
0.08 Myr-i for BH-NS, and 15.3-0.4 Myr"! for BH-BH 
systems (submodel A-B). For sub-solar metallicity NS-NS 
(65.9-6.9 Myr-i) and BH-BH systems (56.7-16.1 Myr^^) 
strongly dominate the DCO merger rate, with relatively 
insignificant rates for BH-NS systems (0.5-0.4 Myr~^). 

The most notable aspect of this variation is the highest 
Galactic NS-NS merger rate: 65.9 Myr^^ for submodel A 
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and sub-solar nictallicity. When compared to the stan- 
dard model, these systems merge ~ 8 times more often 
(from 8.1 Myr~^ to 65.9 Myr~^). In this variation the A 
is smaller (by at least factor of 2) than the one used in the 
standard model for a NS progenitor (typically a 10 M0 
star; see Fig. [T]). It means that the envelope is strongly 
bound in V2. As a consequence the binaries will end the 
CE with smaller orbital separations than they would in 
the standard model. This effect is additionally enhanced 
since NS-NS progenitors, for Z = 0.1 Z0, often experi- 
ence two CE events (see NSNSOl in Tab. [S] in standard 
model; a similar channel is also found in V2). A decreased 
semi-major axis increases the chances of survival through 
a SN. The vast majority of disruptions occur correspond- 
ing to the first SN, when the systems are still wide. For 
example, in the standard model ^ 94% of disruptions are 
encountered at the first SN, while ^ 6% at the second 
SN. In V2, there are virtually no disruptions at the sec- 
ond SN leading to the 8-fold increase of merger rates. A 
secondary effect of the strongly bound envelope is the mi- 
gration of the progenitors toward low merger times. This 
has two competing consequences. One is the increased 
number of mergers during CE and a moderate reduction 
of merger rates. The second is the efficient dissipation of 
orbital energy that allows wide progenitors (that in the 
standard model would produce non-merging NS-NS) to 
form merging NS-NS systems moderately increasing the 
merger rates. 

For solar metallicity the effect of SN survival does not 
manifest itself because progenitors of NSs often experi- 
ence ECSe (high wind mass loss reduces prc-SN progen- 
itor mass), which generally do not disrupt binaries (no 
kick). However, for there is a significant NS-NS merger 
rate drop in submodel B (from 7.6 Myr~^ in the standard 
model to 1.1 Myr~^ in this variation). The main forma- 
tion channel in the standard model (see NSNSOl in Table 
S]) involves a CE and non-conservative mass transfer after 
the first SN. In V2 after the CE the orbital separation is 
smaller than in the standard model. Therefore the sec- 
ondary, a low mass helium star that is a progenitor of the 
second NS, experiences RLOF much earlier in its evolu- 
tion in V2. A combination of shorter orbital period and 
the earlier evolutionary stage (typically helium HG) leads 
to the development of the second CE in V2 (as opposed to 
stable mass transfer in the standard model) . Since the sur- 
vival through CE with HG donors is not allowed in model 
B, the merger rates significantly decrease. There is also a 
small decrease of merger rates in submodel A. Although 
HG CE is allowed in submodel A, small orbital separation 
in V2 sometimes leads to a merger in the second CE event. 

The most notable change in the merger rates of BH- 
NS systems is the drop for Z = 0.1 Zq: 3.4-0.5 Myr'^ 
(standard-this variation) for submodel A and 2.3-0.4 
Myr^^ for submodel B. This comes from the fact that, 
for low metallicity environments, the progenitors of BH- 
NS systems quite often engage in a CE at the beginning 
of their evolution (see BHNSOl in Tab. [5]). In this early 
CE event, the donor is the more massive component and 
a progenitor of the BH. The A value used in this varia- 
tion (0.1) is higher than in the standard model for a BH 
progenitor (see Fig. I3l4p . This means that the envelope 
is easier to eject and after the first CE the binaries will 



not lose as much orbital energy remaining above the 10 
Gyr merger time. Additionally a wider separation after 
the first CE increases the chances of disruption by the 
supcrnovae. A combination of both effects accounts for 
the drop in BH-NS merger rates for both submodels for 
Z = 0.1 Zq. This does not manifest itself for the Z = Zq 
model. High wind mass loss rates present at solar metal- 
licity often reduce the initial mass ratio of the progenitors 
(> 5-^) and the first interaction instead of being a CE is 
a non-conservative mass transfer (see BHNSOl in Tab. |4]). 
Since A plays role on in CE events, the BH-NS rates are 
comparable in V2 and the standard model. The factor of 
~ 2 difference arises from the later CE events initiated by 
the NS progenitor. 

For solar metallicity the BH-BH merger rates increase 
by a factor of 2 (from 8.2 Myr"^ to 15.3 Myr"^) for sub- 
model A, and decrease by factor of ^ 5 (from 1.9 Myr""'^ to 
0.4 Myr^^) for submodel B, in contrast with the standard 
model. The A value fixed in this variation (at 0.1) is larger 
than the average Nanjing values for typical BH progenitors 
(see Figs. [3] and H]). High A and correspondingly low bind- 
ing energy leads to small orbital contraction during the CE 
phase. This can (counter-intuitively) both increase and 
decrease merger rates of various BH-BH populations. For 
example, consider BH-BH populations in submodel A and 
B. From the relative merger rates (A much higher than B) 
it is clear that submodel A is dominated by BH-BH bina- 
ries that evolved through the CE phase with a HG donor. 
In contrast, submodel B includes only systems that have 
formed through CE with a CHeB donor (by assumption). 
Since HG stars are smaller (for the same mass and metal- 
licity) than CHeB stars, progenitors in submodel A have 
relatively small orbital separations, while in submodel B 
they have relatively wide separations. On one hand, short 
period progenitor binaries do not merge in the CE phase 
as often as in the standard model, since the A chosen in 
V2 is large, and that leads to an increased BH-BH merger 
rate in submodel A. On the other hand, long period pro- 
genitor binaries are not sufficiently contracted (in term of 
separation) due to the high A in V2. They form long- 
period BH-BH systems with merger times exceeding 10 
Gyr, and thus the merger rate decreases in submodel B. 
For sub-solar metallicity the effects of an increased A are 
compensated for by the smaller radial expansion of stars as 
compared to solar metallicity, and the V2 merger rates for 
both submodels are similar to the standard model rates. 

The distribution of chirp masses for both submodels and 
metallicities closely resembles that of the standard model. 
The average chirp masses for Z = Zq arc 1.08-1.06 Mq 
(submodel A-B) for NS-NS, 3.2-3.2 Mq for BH-NS, and 
6.5-6.5 Mq for BH-BH systems. For Z 0.1 Zq the cor- 
responding values are 1.09-1.07 Mq (submodel A-B) for 
NS-NS, 3.5-3.3 Mq for BH-NS, and 17.2-9.3 Mq for BH- 
BH systems. 

The average values of delay times for Zq are: 90-618 
Myr (submodel A-B) for NS-NS, 2.2-2.6 Gyr for BH-NS, 
and 1.8-6.5 Gyr for BH-BH systems. For 0.1 Zq the values 
are: 0.2-1.6 Gyr (NS-NS), 1.6-1.7 Gyr (BH-NS), 1.2-1.9 
Gyr (BH-BH systems). 



4.4. Variation 3 
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This variation fixes the binding to A = 1. The Galac- 
tic merger rates for are dominated by NS-NS systems 
(48.8-14.3 Myr-\ submodel A~B) followed by BH-BH 
(5.0-0.03 Myr-i ) and BH-NS systems (4.6-0.03 Myr"!). 
For sub-solar metallicity, the BH-BH systems merge most 
often (90.2-7.9 Myr-\ submodel A-B), then NS-NS sys- 
tems (44.1-4.2 Myr-i) and BH-NS systems (15.8-8.4 
Myr-i ). 

For solar metallicity there are no significant changes in 
merger rates for submodel A. However, the rise in rates for 
NS-NS systems by a factor of ^ 2 when compared to the 
standard model (from 23.5 Myr~^ to 48.8 Myr~^) makes 
it the most frequently merging population of double NSs, 
for Zq, in this study. The increase of merger rates for NS- 
NS and BH-NS systems (from 1.6 Myr'^ to 4.6 Myr^^) is 
associated with fewer binaries merging during the CE due 
to the high A value used in this variation. However, the 
high A values have the opposite effect for BH-BH systems. 
The resulting low binding energy of the CE is responsi- 
ble for the slight drop in merger rates for BH-BH systems 
(from 8.2 Myr~^ to 5.0 Myr"-'^), as fewer binaries cross the 
10 Gyr point (toward the merging population) due to in- 
sufficient orbital energy loss when ejecting their envelopes. 

The population of NS-NS systems in submodel B follows 
the behaviour described in submodel A. What is noticeable 
for BH-NS and BH-BH systems in submodel B. at Zq, is 
a drop in the merger rates (by an order of magnitude with 
respect to the standard model) . By assumption, submodel 
B favours DGOs which initiated a CE with a donor beyond 
the HG, in this case during CHeB. Donors of this type that 
are progenitors of components in BH-NS or BH-BH sys- 
tems have lower (see Figs. [2]-[4|) A values in the standard 
model (more physical) than in this variation (A = 1). A di- 
rect consequence is the reduced orbital energy loss during 
a CE event. This allows binaries with CHeB CE to retain 
a large separation at the end of the CE phase, preventing a 
DCO merger within 10 Gyr. Therefore, the population of 
merging DCOs containing BHs for submodel B is reduced. 

For sub-solar metallicities the high A effect is still rel- 
evant. However, in these chemical environments donors 
initiate CE more often as CHeB stars (instead of HG) due 
to smaller radial expansion during the evolution. This in- 
creases the number and merger rates of NS-NS and BH-NS 
systems in submodel B. This also allows for more BH-BH 
systems to be formed in this submodel but for this type 
of DCO the high A causes a counter-effect, and efficiently 
reduces the merger rates of double BHs. 

The average chirp masses for Z0 arc: 1.06-1.05 Mq 
(submodel A-B) for NS-NS, 2.7-2.4 Mq for BH-NS, and 
6.0-5.9 Mq for BH-BH systems. For 0.1 Zq the values 
are: 1.09-1.12 Mq for NS-NS, 2.9-2.9 Mq for BH-NS, and 
12.5-6.8 for BH-BH systems. 

The average delay times for Zq are: 1.2-2.2 Gyr (sub- 
model A-B) for NS-NS, 2.0-2.7 Gyr for BH-NS, and 4.2- 
2.2 Gyr for BH-BH systems. The values for 0.1 Zq are: 
0.9-2.4 Gyr for NS-NS, 1.8-2.3 Gyr for BH-NS, and 1.6-3.3 
Gyr for BH-BH systems. These average values are greater 
than for the standard model, as the high A value used in 
this variation (low binding energy) allows the binaries to 
retain a larger orbital separation after the CE phase. 

4.5. Variation 4 



In this model we fix the A value at 10. The Galactic 
merger rates for Zq are: 20.8-0.3 Myr~^ (submodel A-B) 
for NS-NS, 0.9-0.0 Myr"! for BH-NS, and 0.3-0.0 Myr"! 
for BH-BH systems. The rates for 0.1 Zq are: 29.5-1.4 
Myr-i for NS-NS, 8.8-1.6 Myr^i for BH-NS, and 5.9-0.3 
Myr"i for BH-BH systems. 

For Zq the populations of DCOs of all types experi- 
ence a reduction in the merger rate for both submodels 
(see Table [2]) when compared with the standard model. 
The mechanism responsible for this result is straightfor- 
ward. The very high A value means a low binding energy 
of a star's envelope, and so it is easily ejected during a 
CE event. As a consequence, little orbital energy is lost, 
and the separation of the components is reduced insignif- 
icantly. The weak orbital contraction causes the DCOs 
to retain merger times larger than 10 Gyr, thus creating 
a large non-merging population at the cost of the merg- 
ing one. The mechanism described above holds true for 
all types of DCOs, but manifests most clearly for BH-BH 
systems. For example, for submodel A, for each BH-BH 
system formed below the 10 Gyr merger time limit, ~ 240 
are formed above it. For submodel B the rates are null for 
BH-NS and BH-BH systems. In other words, despite high 
statistics — 6 million binary stars simulated — no DCO sys- 
tems containing BHs were produced. In order to under- 
stand this result one needs to recall that the CE phase 
is usually initiated by HG and CHeB donors, due to sig- 
nificant radial expansion during these stages of evolution. 
Submodel B does not allow for a HG donor so the only 
systems left to populate this submodel are those with a 
CHcB donor during the CE, in their evolutionary history. 
This evolved CE donor indicates that these binaries had 
initial separations large enough to be able to bypass the 
CE on the HG, and this initial orbital condition places 
them in the long merger time regime. For progenitors of 
BH-NS and BH-BH systems the A value in the standard 
model during CHeB is much lower (see Figs. [2H11) when 
compared to the one used in this variation (10). So, due to 
insufficient orbital energy loss during CE, no binaries are 
drawn toward the merging population, making it void of 
system containing BHs. This migration of DCOs between 
the merging and non-merging populations is illustrated is 
Figure [HI 

For 0.1 Zq more merging BH-NS/BH systems arc 
formed in general due to weaker SNe (see Sec. 14. 1|) . The 
probability of producing such a system for the merging 
population increases, yielding non-zero merger rates for 
submodel B. Additionally, the merger rates for BH-NS sys- 
tem for submodel A increase by a factor of ~ 3 (from 3.4 
to 8.8 Myr^^) when contrasted with the standard model. 
This result is opposite to the one for Zq , where the merger 
rate of BH-NS systems for submodel A drops slightly (from 
1.6 to 0.9 Myr-i ). A weakly bound CE (in this varia- 
tion) causes insignificant orbital reduction and increases 
the chances of disruption by the following SNe (after the 
first or second CE, sec Table 0. However the stalled SN 
kicks associated with low metallicity environments reduce 
the disruption probability, accounting for the increase in 
0.1 Zq in V4. 

For NS-NS systems the Galactic merger rates increase 
from 8.1 Myr"^ (standard model) to 29.5 Myr"^ (this 
variation) at 0.1 Zq. A characteristic of these DCOs for 
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sub-solar metallicitics is the second CE initiated just after 
the first one (see Table . Two weakly bound CE phases 
allow the binary to avoid a merger in either of these events, 
effectively increasing the number and merger rates of NS- 
NS systems. 

The average chirp masses for Zq are ~ 1.03 M0 for both 
submodels for NS-NS systems. For BH-NS and BH-BH 
systems, submodel A, these are: 2.5 Mq and 5.8 Mq, re- 
spectively. For 0.1 Zq these values are: 1.1-1.2 M0 (sub- 
model A-B) for NS-NS, 2.9-2.7 Mq for BH-NS, and 7.6- 
6.7 Mq for BH-BH systems. 

The average delay time values for Zq, for NS-NS sys- 
tems are 0.5 Gyr and 2.0 Gyr, for submodel A and B, 
respectively. For BH-NS and BH-BH systems these are 
2.5 Gyr and 3.2 Gyr, respectively, for submodel A. For 
0.1 Zq these values are: 0.7-2.0 Gyr (submodel A-B) for 
NS-NS, 2.4-2.2 Gyr for BH-NS, and 3.8-2.7 for BH-BH 
systems. 

4.6. Variation 5 & 6 

In variations 5 and 6 we set the maximum mass achiev- 
able for a NS to 3.0 Mq and 2.0 Mq, respectively (this 
is 2.5 Mq in the standard model). These models produce 
no statistically significant changes, when compared to the 
standard model. The reason for the population of DCOs 
being insensitive to the maximum mass of a NS is the usage 
of the Rapid supernova mechanism. The main character- 
istic of this engine is the lack of compact objects formed 
in the range ~ 2-5 Mq. 

The distribution of delay times and chirp masses, as well 
as their average values, are nearly identical (within statis- 
tical accuracy) to the standard model. 

4.7. Variation 7 

In this variation we set the a of the Maxwellian distri- 
bution of the natal kicks to 132.5 km/s (it is 265 km/s in 
the standard model). The Galactic merger rates for Zq 
are: 32.4-9.5 Myr"! (submodel A-B) for NS-NS, 1.9- 
0.3 Myr-i for BH-NS, and 10.4-2.1 Myr"! for BH-BH 
systems. For sub-solar metallicity the rates are: 8.3-2.2 
Myr-i for NS-NS, 6.1-4.3 Myr"! for BH-NS, and 83.7- 
15.1 Myr-i . 

For solar metallicity we find slightly increased (by less 
than a factor of ^ 2) Galactic rates for all DCO types, 
for both submodels. In practice, reducing the a value 
means that natal kicks of high magnitude are chosen less 
frequently from the Maxwellian distribution. As a conse- 
quence fewer binaries are disrupted during the SN explo- 
sions. This increases the pool of potential merging DCOs, 
and accounts for higher merger rates in this variation. The 
least affected, when compared to the standard model, are 
BH-BH systems, since the progenitors of these systems are 
massive stars, and thus produce weak SNe (stalled or no 
kicks, direct collapse). Therefore, statistically, these sys- 
tems are minimally impacted by decreasing the occurrence 
of high magnitude natal kicks in the binary population. 

For sub-solar metallicity this effect is even less signifi- 
cant, due to naturally weaker explosions. Low wind mass 
loss rates associated with low metallicity environments 
produce more massive pre-SN progenitors. These stars 
cause more fallback of ejecta during the explosion, due to 
their strong gravitational potential. Therefore, the mag- 
nitude of the velocity of the kicks is in general smaller 



than for Zq . Further reduction of the occurrence of high 
magnitude kicks yields no statistically significant effects. 

In the case of NS-NS systems for 0.1 Zq the rates are al- 
most unchanged when compared with the standard model. 
Note that most of these DCOs have two CE events (ex- 
actly like in the standard model, see Table[5]), which makes 
them susceptible to merging during one of these episodes. 
Reducing the probability of high velocity kicks results in 
smaller separations after the first SN. A consequence is an 
increased number of CE mergers that counters increased 
SN survivability. 

The average chirp masses for Zq are: 1.05 Mq for NS- 
NS, 3.1 Mq for BH-NS, and 6.6 for BH-BH systems (aU for 
both submodels). For 0.1 Zq the values are: 1.08 Mq for 
NS-NS and 3.1 Mq for BH-NS systems (for both submod- 
els). For BH-BH systems these are: 12.4 Mq for submodel 
A and 9.7 Mq for submodel B. 

The average values of delay time distribution are: 1.3- 
2.2 Gyr (submodel A-B) for NS-NS, 2.7-4.3 Gyr for BH- 
NS, and 1.2-2.9 Gyr for BH-BH systems. For 0.1 Zq the 
values are: 1.0-2.7 Gyr for NS-NS, 1.9-1.9 Gyr for BH-NS 
and 1.1-1.3 Gyr for BH-BH systems. 

4.8. Variation 8 

In this model we only enforce full natal kicks for the 
BHs. The results therefore only change for BH-BH/NS 
systems. The Galactic merger rates for Zq for BH-NS 
systems arc 0.03-0.004 Myr'^ (submodel A-B). For BH- 
BH systems: 0.05-0.005 Myr'^ For 0.1 Zq the rates are: 
0.7-0.2 Myr-i for BH-NS and 4.2-0.8 Myr"! for BH-BH 
systems. 

The disruption rate of binary progenitors of BH-BH/NS 
systems during the SN, due to the strong velocity kicks, is 
significant. When contrasted with the standard model, the 
Galactic merger rates for Zq drop by 2-3 orders of magni- 
tude (see Table (2). For example, for BH-BH, submodel B 
at Zq, the merger rate drops from 1.9 Myr~^ (standard) 
to 0.005 Myr~\ For 0.1 Z© the effect is smaller but stiU 
noticeable, as the rates drop by a factor of ^ 10 (see Ta- 
ble[3]). The pre-SN progenitors of BHs in this environment 
are, on average, more massive when compared to Zq. This 
results in strongly bound binaries, which makes them less 
susceptible to being disrupted from natal kicks. 

The average chirp masses for Zq are: 3.0-3.2 Mq (sub- 
model A-B) for BH-NS systems. For BH-BH systems 
these are 6.5 Mq and 5.9 Mq for submodels A and B, 
respectively. For 0.1 Zq these values are 1.09 Mq (both 
submodels) for NS-NS systems and 3.0 Mq (also both sub- 
models) for BH-NS systems. For BH-BH systems these are 
9.0 Mq and 7.2 Mq for submodels A and B, respectively. 

The mean delay times for Zq are: 1.2-1.7 Gyr (sub- 
model A-B) for NS-NS systems and 2.0-7.0 Gyr for BH- 
BH systems. For BH-NS, submodel A, the mean value 
is 0.4 Gyr. For 0.1 Zq the mean values are: 0.9-2.2 Gyr 
(submodel A-B) for NS-NS, 1.2-1.3 Gyr for BH-NS, and 
0.6-1.2 Gyr for BH-BH systems. 

4.9. Variation 9 

In this variation BHs do not receive natal kicks. We 
therefore describe the results for BH-NS and BH-BH sys- 
tems (NS-NS systems behave as in the standard model). 
The Galactic merger rates for BH-NS systems are 1.4- 
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0.2 Myr^i (submodel A-B) for Zg and 5.2-3.7 Myr^^ 
(0.1 Zg). For BH-BH systems the rates are: 16.9-4.2 
Myr-i ( Zg) and 92.3-19.3 Myr"! (0.1 Zg). 

For Zg the merger rates for BH-BH systems mcrease 
by a factor of ^ 2 for both submodels, when compared 
to the standard model. This is a straightforward conse- 
quence of the assumption of no natal kicks used in this 
variation. The survival rate of BH-BH/NS systems in- 
creases due to fewer disruptions occurring during the SN. 
This allows more DCOs with BHs, which in turn increases 
the merger rates. For 0.1 Zg the gain is smaller. Pre-SN 
progenitors of BHs are heavier in this chemical environ- 
ment due to low wind mass loss. This causes a lot of 
ejecta fallback (almost all matter is pulled back onto the 
newly formed BH). Setting no natal kicks for these sys- 
tems changes little as these BHs "naturally" suppress the 
SN explosion. 

The behaviour of BH-NS systems is explained analo- 
gously to BH-BH DCOs. However the effects are smaller 
since the survival rates due to disruption are dominated 
by SNe resulting in the formation of NSs. Pre-SN progeni- 
tors of NSs arc light and experience low fallback (stronger 
kicks). This allows for more powerful explosions, which 
disrupt the binaries and prevent the creation of DCOs. 

The average chirp masses for BH-NS systems are: 
3.1 Mg for Zg and 3.0 Mg for 0.1 Zg (all for both sub- 
models). For BH-BH systems the values are: 6.3-6.2 Mg 
(submodel A-B) for Zg and 12.1-9.3 Mg for 0.1 Zg. 

The mean values of delay times for BH-NS systems 
are: 1.7-3.0 Gyr (submodel A-B) for Zg and 1.7-1.7 Gyr 
(0.1 Zg). For BH-BH systems these are: 1.5-3.4 Gyr 
for Zg and 1.2-1.4 Gyr for 0.1 Zg 

4.10. Variation 10 

In this variation we use the Delayed supernova engine 
(as compared with the Rapid engine used in the standard 
model). The Galactic merger rates for Zg are: 25.6-8.9 
Myr-i (submodel A-B) for NS-NS, 0.07-0.03 Myr"! for 
BH-NS, and 0.6-0.08 Myr"! for BH-BH systems. For 
0.1 Zg the rates are: 8.6-2.6 Myr"! for NS-NS, 2.3-2.0 
Myr-i for BH-NS, and 62.0-11.5 Myr"! for BH-BH sys- 
tems. 

For all metallicities and submodels we note a drop in 
the merger rates for systems containing BHs (see Table 
ID and [3]). The reason for this is the difference in the 
amount of fallback matter acquired in SN explosions re- 
sulting from the two engines. To illustrate, we present 
the resu lts of an example c alculation using Eq.l6, 18, and 
19 from lFrver et al.l (|2011h . estimating the fallback factor 
/tb. For a star prior to the SN explosion with a total mass 
M = 8.7 Mg and carbon-oxygen core Mqo = 6.5 Mg 
(typical core mass collapsing to a BH), we get the fol- 
lowing values of fpo' 1-0 (full fallback) for the Rapid and 
0.46 (partial fallback) for the Delayed engine. For progen- 
itors of BHs the amount of fallback is lower (more mass 
is ejected) in the case of the Delayed engine. This results 
in natal kicks of a higher magnitude (see Eq. [3]), which 
are more capable of disrupting the binary and preventing 
the creation of a DCO. For progenitors of NSs the fall- 
back is similar in both VIO and the standard model, and 
therefore we find similar magnitudes of natal kicks, and a 
similar chance for disruption by SN. 

Another significant result produced by this model is the 



lower limit of the chirp mass distribution for BH-BH sys- 
tems: ^ 2 (for an example see Figure fTU]) . The Delayed 
supernova engine creates remnant objects (NSs and BHs) 
with a continuous spectrum of masses, contra r y to t he 
Rapid engine (for details see iBelczvnski et al.l ( 2011b[ )). 
Thi s suggests observational tests to distinguish the mod- 
els (jFrver et al.U201lUBelczvnski et all , feoilbl) . 

The average chirp mass for Zg is: 1.14-1.13 Mg (sub- 
model A-B) for NS-NS, 3.1-3.0 Mg for BH-NS, and 5.7- 
4.6 Mg for BH-BH systems. For 0.1 Zg these are: 1.2- 
1.2 Mg for NS-NS, 3.4-3.4 Mg for BH-NS, and 14.4- 
10.2 Mg for BH-BH systems. 

The mean delay time values for Zg are: 1.2-1.8 Gyr 
(submodel A-B) for NS-NS, 2.5-1.4 Gyr for BH-NS, and 

1.7- 2.6 Gyr for BH-BH systems. For 0.1 Zg these are: 
0.9-2.1 Gyr for NS-NS, 1.8-1.7 Gyr for BH-NS, and 1.2- 
1.5 Gyr for BH-BH systems. 

4.11. Variation 11 

In this variation we reduce the stellar wind mass loss 
rates by a factor of 2. The Galactic merger rates for Zg 
are: 24.2-6.5 Myr"! (submodel A-B) for NS-NS, 1.2-0.2 
Myr-i for BH-NS, and 29.7-3.6 Myr"! for BH-BH sys- 
tems. For 0.1 Zg the rates are: 7.7-2.3 Myr"! for NS-NS, 

3.8- 2.4 Myr-i for BH-NS, and 79.2-17.1 for BH-BH sys- 
tems. 

The most significant result due to this variation is the 
increase in the Galactic merger rate for BH-BH systems 
for Zg. The rates increase by a factor of ~ 3 and ~ 2 
for submodels A and B, respectively. The effect of the 
lowered wind mass loss rates on double BHs is manifold. 
First, the decreased reduction of mass of BH progenitors 
causes them to retain more mass just prior to the SN. 
This results in a larger amount of fallback matter during 
the explosion, which stalls the natal kicks. This in turn 
increases the survivability rate of binaries and produces 
more merging systems. Second, the increased prc-SN mass 
also results in more massive remnants, which is clearly vis- 
ible in Figure [TU] and 1111 Third, the minimum mass of a 
BH progenitor is reduced to lower values (from the usual 
~ 20 Mg to ^ 18 Mg). For sub-solar environments these 
effects are weak as the wind mass loss rate for 0.1 Zg is 
already low. Decreasing this even further yields no new 
qualitative results. 

For NS-NS systems these effects are also insignificant. 
The progenitors of NSs in these systems are lighter than 
those of BHs, and for these stars the wind mass loss rates 
are smaller than for their heavier counterparts. Therefore, 
decreasing wind mass loss even further plays no significant 
role in the survivability or the properties of these systems. 

The merger rates for BH-NS systems do not increase in 
the way they do for BH-BH systems, despite their har- 
boring a BH and a massive NS progenitor. Binary stars 
that would have formed BH-NS systems in the case of 
"standard" winds, now retain more mass during evolution 
(due to the reduced winds), and form BH-BH systems in- 
stead. This reduces the BH-NS population and adds merg- 
ing DCOs to the BH-BH systems. 

The average chirp masses for Zg are: 1.05-1.05 Mg 
(submodel A-B) for NS-NS, 3.2-3.3 Mg for BH-NS, and 
10.5-9.1 Mg for BH-BH systems. Note that for the latter 
DCOs the range of chirp mass is the largest (~ 14 Mg) for 
this metallicity of all variations. For 0.1 Zg the average 
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chirp masses arc: 1.08-1.08 Mq for NS-NS, 3.3-3.4 Mq 
for BH-NS, and 14.3-10.6 Mq for BH-BH systems. 

The mean delay times for Zq are: 1.0-1.6 Gyr (sub- 
model A-B) for NS-NS, 1.3-1.9 Gyr for BH-NS, and 0.6- 
2.8 Gyr for BH-BH systems. For 0.1 Z© these values are: 
1.0-2.3 Gyr for NS-NS, 1.5-1.5 Gyr for BH-NS, and 1.0- 
1.2 Gyr for BH-BH systems. 

5. SUMMARY 

We have investigated the major parameters and input 
physics involved in binary evolution leading to the for- 
mation of double compact objects: the CE coefficient A, 
the supernova mechanisms (remnant formation and na- 
tal kick magnitude), the maximum mass of NSs, and the 
wind mass loss rates. The study was performed by calcu- 
lating a suite of population synthesis models, allowing us 
to estimate the associated uncertainties in the formation 
of DCOs. 

The calculations were done using the Star Track pop- 
ulation synthesis code, recently updated with wind mass 
loss rates for massive stars, a realistic CE treatment, and 
convection-driven neutrino-enhanced supernovae. These 
updates are incorporated in our revised Standard Model 
(see Section 14. ip of binary evolution and DCO forma- 
tion. Our newer version yields lower merger rates than the 
equivalent Standard Model from .Bclc zynski et al. (200_2l). 
This difference arises mainly from the Nanji ng treat- 
ment of the CE parameter, A (see Section 12.3.21 IXu fc Lil 
(|2010D 1. which we have now incorporated. Significant 
changes in the merger rates are also caused by varying 
the value of A (Variations 1-4, see Table [2] and [3]) . We 
therefore identify a strong dependence of the merger rates 
on the binding energy of the CE. 

Emerging from the simulations are two populations of 
DCOs. The first is that of double compact objects with 
merger times less than 10 Gyr (merging in a Milky Way- 
like galaxy, the merging population). The second is of 
objects having merger times higher than this limit (dis- 
tant future mergers, the non-merging population). The 
A parameter shifts DCOs between these two populations. 
On the one hand, lowering A increases the CE binding 
energy, and thus tends to draw the binaries from the non- 
merging population toward the merging one, as progeni- 
tor binaries must lose more orbital energy to unbind the 
envelope. This results in a larger number of system in 
the merging population, and increases the Galactic merger 
rates. On the other hand, this also causes progenitors of 
DCOs from the merging population to merge during the 
CE more often, terminating binary evolution and therefore 
reducing the number of binary systems, and thus merger 
rates. The final outcome is set by the specific value of A, 
and varies from population to population. For example, 
a given change in A may cause an increase in the merger 
rate of BH-BH systems while decreasing the rate of NS-NS 
systems, due to the differing orbital separations of BH-BH 
and NS-NS progenitors. 

The types of DCOs most affected by this mechanism are 
NS-NS and BH-NS systems. The range of merger rates for 
these systems spans three order of magnitude between so- 
lar and sub-solar metallicities (sec Table [2] and Table [3]). 
The large range of merger rates obtained in this and pre- 
vious studies using calculations with a fixed value of A is 



an artifact of the unphysical treatment of the CE phase. 
When a more physically motivated treatment is applied, 
in which the value of A is chosen to correspond to the given 
star and its evolutionary stage, it is found that the merger 
rates are on the high side of the ranges computed with 
variations that assume non-physical, fixed values of A 

In particular, for the standard model, for either choice 
of metallicity of submodel (A or B), the BH-BH Galactic 
merger rate is at least 1.9 per Myr. This is a factor of ~ 5 
higher than the "realistic" estimate used in lAbadie et al.l 
(j201Clt) , which translates into a detection rate of 20 BH-BH 
mergers per year with Advanced LIGO and Virgo detec- 
tors in that paper. Therefore, if the standard model is 
correct, we may expect that advanced gravitational-wave 
detectors will be able to capture above a hundred BH- 
BH coalescences per year. This is a clearly positive pre- 
diction from the perspective of the ongoing searches for 
gravitational-wave signals. However, we are unable to pro- 
vide error estimate on this improvement, since there are 
no alternative physical estimates for A readily available to 
us. Therefore, one should treat this estimate with some 
caution. 

We can test the validity of certain models using 
the observed limits on the Galactic NS-NS merger 
rate ()Kim et al.l . |2006[) . Their upper limit (190 Myr"!) 
is not violated by any of our models. However, the lower 
limit (3 Myr~^) is not met by several of our models, and 
they would therefore be ruled out. Models (for Zq) that do 
not meet this condition are: Variation 1-submodel A (NS- 
NS Galactic merger rate 0.4 Myr~^), Variation 1-submodel 
B (0.4 Myr-i), Variation 2-submodel B (1.1 Myr'^), and 
Variation 4-submodel B (0.3 Myr^^); this is illustrated in 
Figure [121 All these variations employ a fixed value of A, 
and the observations of Galactic NS-NS systems lends fur- 
ther support to our claim that such a treatment is unphysi- 
cal. Additionally we can use the estimates on the predicted 
chirp mass of ICIO X-1 (0.3 Zq environment) to further 
validate our results. We test our chirp mass distribution 
for 0.1 Zq against the predicted mini mal predicted chri p 
mass for ICIO X-1 equal to 15 Mq (jBulik et al.l . |2011[) . 
From Fig. [13] and Table [9| one can see that models V3- 
submodcl B (max. chirp mass 14.9 Mq) and V4-submodel 
B (max. chirp mass 11.5 Mq) do not reach this limit. 
Again these models use constant, in this case high A val- 
ues. This is further evidence to support the claim that 
such a treatment is unpshysical. 

Another infiuential factor in DCO formation is the treat- 
ment of SN explosions. Specifically, the magnitude of the 
natal kick plays a significant role, as the high kicks (as- 
sumed in our standard model) tend to disrupt DCO pro- 
genitors, reducing the merger rates. The rates for reduced 
(by half) kicks are given in model V7. The correspond- 
ing increase in the rates is rather modest, due to the fact 
that for NS-NS systems some NSs explode as electron cap- 
ture supernovae and do not receive kicks at all. BH kicks 
are already reduced by the amount of fall back and further 
kick decrease (coming from reducing a in V7) does not sig- 
nificantly increase the rates. However, the effect of kicks 
becomes significant for BH-BH systems once the assump- 
tion on the reduction of kicks by fallback is relaxed. We 
have allowed for the very wide range of possible BH kicks 
in models V8 (high kicks) and V9 (no kicks). It is found 
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that the Galactic merger rates of BH-BH systems change 
by ^ 2-3 orders of magnitude (see Table [51 [3]) . Therefore, 
the mostly unknown magnitude of BH natal kick limits 
our predictive power for the BH-BH rate estimates. 

Despite the fact that we lack strong observational or 
theoretical constraints on BH natal kicks, in our standard 
model we have adopted the most likely model of natal 
kicks (ones that decrease with increasing BH mass). This 
model is supported by existing observations (most mas- 
sive BHs seem to form without a kick) and can be intu- 
itively understood in the framework of natal kicks arising 
from asymmetric supernova mass ejection. The standard 
model BH-BH merger rates are close to upper limit set by 
the full allowed range of possibilities (V8-V9), and again 
it appears that our more realistic treatment of DCO for- 
mation favors higher rates, and sets an optimistic tone for 
near-future gravitational wave inspiral detection. 

This the first in a series of three pa- 
pers. Results presented here, and their exten- 
sions in the Synthetic Universe on-line database 
(www.syntheticuniverse.org), will be used in the sec- 



ond study, where we will investigate the NS-NS/BH- 
NS/BH-BH merger rates as a function of cosmological 
redshift, star formation rate, and metallicity. The third 
paper will focus on gravitational wave detection rates 
for upcoming observatories (advanced LIGO/Virgo and 
Einstein Telescope). 
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Table 1 
Summary of Models'^ 



Model 


Parameter 


Description 


S 


Standard 


A ^Nanjing, il/Ns,max = 2.5 Mq, cr = 265 
km s^^ BH kicks: variable, SN: Rapid 


VI 


A = 0.01 


very low A: fixed 


V2 


A = 0.1 


low A: fixed 


V3 


A = 1 


high A: fixed 


V4 


A = 10 


very high A: fixed 


V5 


A/NS,max = 3.0 Mq 


high maximum NS mass 


V6 


MNS^max = 2.0 Mq 


low maximum NS mass 


V7 


cr 132.5 km s"! 


low kicks: NS/BH 


V8 


full BH kicks 


high natal kicks: BH 


V9 


no BH kicks 


no natal kicks: BH 


VIO 


Delayed SN 


NS/BH formation: Delayed SN engine 


VII 


weak winds 


Wind mass loss rates reduced to 50% 



All parameters, except for the one listed under "Description", retain 
their Standard model ("S") values. 
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Table 2 

Galactic Merger Rates, Zq [ Myr~^] 



Model 


NS-NS 


BH-NS 


BH-BH 


s 


23.5 (7.6) 


1.6 (0.2) 


8.2 (1.9) 


VI 


0.4 (0.4) 


0.002 (0.002) 


1.1 (1.1) 


V2 


11.8 (1.1) 


2.4 (0.08) 


15.3 (0.4) 


V3 


48.8 (14.3) 


4.6 (0.03) 


5.0 (0.03) 


V4 


20.8 (0.3) 


0.9 (0.0) 


0.3 (0.0) 


V5 


24.1 (8.1) 


1.4 (0.2) 


8.3 (2.0) 


V6 


24.1 (8.3) 


1.4 (0.2) 


8.0 (1.9) 


V7 


32.4 (9.5) 


1.9 (0.3) 


10.4 (2.1) 


V8 


23.3 (7.7) 


0.03 (0.004) 


0.05 (0.005) 


V9 


23.4 (8.0) 


1.4 (0.2) 


16.9 (4.2) 


VIO 


25.6 (8.9) 


0.07 (0.03) 


0.6 (0.08) 


VI 1 


24.2 (6.5) 


1.2 (0.2) 


29.7 (3.6) 


Range 


0.4-48.8 (0.4-9.5) 


0.002-4.6 (0.0-0.3) 


0.05-29.7 (0.0-4.2) 



^ Rates are calculated for a synthetic galaxy similar to the Milky Way 
(solar metallicity, and 10 Gyr of continuous star formation at the level of 
3.5 M0 yr^^). Rates are presented for submodel A (CE HG donor allowed), 
with the rates for submodel B listed in parentheses (CE HG donor forbidden); 
see Section 12.3.11 for details. The range presents the minimum and maximum 
value for each DCO type. 



Table 3 

Galactic Merger Rates, 0.1 Z© [ Myr""^] 





NS-NS 


BH-NS 


BH-BH 


s 


8.1 (2.5) 


3.4 (2.3) 


73.3 (13.6) 


VI 


0.06 (0.06) 


0.03 (0.03) 


12.5 (8.1) 


V2 


65.9 (6.9) 


0.5 (0.4) 


56.7 (16.1) 


V3 


44.1 (4.2) 


15.8 (8.4) 


90.2 (7.9) 


V4 


29.5 (1.4) 


8.8 (1.6) 


5.9 (0.3) 


V5 


8.0 (2.3) 


3.4 (2.1) 


73.4 (13.7) 


V6 


7.8 (2.4) 


3.5 (2.0) 


74.5 (13.8) 


V7 


8.3 (2.2) 


6.1 (4.3) 


83.7 (15.1) 


V8 


8.2 (2.5) 


0.7 (0.2) 


4.2 (0.8) 


V9 


8.1 (2.1) 


5.2 (3.7) 


92.3 (19.3) 


VIO 


8.6 (2.6) 


2.3 (2.0) 


62.0 (11.5) 


Vll 


7.7 (2.3) 


3.8 (2.4) 


79.2 (17.1) 


Range 


0.06-65.9 (0.06-6.9) 


0.03-15.8 (0.03-8.4) 


4.2-92.3 (0.3-19.3) 



Same as Table [3] but for sub-solar metallicity. 
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Table 4 

Formation channels of DCOs for Zq 



Channel Fraction 



NSNSOl NC:a->b, SN:a, CE:b^a, NC:b^a, SN:b 79.3% 

NSNS02 NC:a->b, CE:b->a, NC:b->a, AIC(WD^NS):a, NC:b->a, SN:b 8.0% 

NSNS03 NC:a^b, SN:a, CE:b^a, CE:b^a^SN:b 6.9% 

NSNS04 Other 5.8% 

BHNSOl NC:a-^b, SN:a, CE:b->a, SN:b 95.4% 

BHNS02 NC:a->b, SN:a, CE:b^a, NC:b^a, SN:b 1.8% 

BHNS03 Other 2.8% 

BHBHOl NC:a^b, SN:a, CE:b^a, SN:b 98.9% 

BHBH02 Other 1.1% 



Coalescing DCOs' formation channels for the standard model, submodel A at solar metallic- 
ity. NC: non-conservative mass transfer, SN: supernova, CE: common envelope, AIC: accretion 
induced collapse of oxygen/neon white dwarf into NS. The arrows show the direction of transfer 
and "a" stands for the primary (initially more massive) component, "b" for the secondary. 

^The first CE is initiated by the H-rich Hertzsprung gap donor (allowed in model A). The 
second starts when the exposed core of the donor becomes an evolved helium star. 



Table 5 

Formation channels of DCOs for 0.1 







Channel 


Fraction 


NSNSOl 


NC 


:a^b, SN:a, CE:b^a, CE:b^a, SN:b 


49.1% 


NSNS02 


NC: 


:a->b, SN:a, CE:b^a, NC:b^a, SN:b 


21.2% 


NSNS03 




NC:a-^b, SN:a, CE:b^a, SN:b 


18.2% 


NSNS04 


NC 


:a^b, CE:b^a, SN:a, CE:b->a, SN:b 


3.3% 


NSNS05 




Other 


8.2% 


BHNSOl 


CE: 


:a^b, SN:a, CE:b->a, NC:b-^a, SN:b 


40.8% 


BHNS02 




CE:a^b, SN:a, CE:b^a, SN:b 


17.4% 


BHNS03 




NC:a->b, SN:a, CE:b^a, SN:b 


13.4% 


BHNS04 


NC: 


:a^b, SN:a, CEib^a, NC:b^a, SN:b 


12.2% 


BHNS05 


NC: 


:a-^b, CE:b^a, NC:a^b, SN:a, SN:b 


8.8% 


BHNS06 




Other 


6.4% 


BHBHOl 




NC:a^b, SN:a, CE:b->a, SN:b 


90.6% 


BHBH02 




CE:a^b, SN:a, CE:b->a, SN:b 


4.0% 


BHBH03 


NC: 


:a->b, SN:a, NC:b^a, CE:b^a, SN:b 


1.4% 


BHBH04 




Other 


4.0% 



^ Same as Table |4] but for sub-solar metallicity. 
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Table 6 

Chirp mass characteristics for Zq, submodel A' 





min 


NS-NS 
avg 


max 


min 


BH-NS 
avg 


max 


min 


BH-BH 

avg 


max 


s 


0.96 


1.05 


1.40 


2.3 


3.2 


3.7 


5.1 


6.7 


8.7 


VI 


1.08 


1.09 


1.14 


3.2 


3.2 


3.2 


5.3 


6.5 


8.3 


V2 


0.96 


1.08 


1.53 


2.5 


3.2 


4.0 


5.0 


6.6 


8.4 


V3 


0.94 


1.06 


1.69 


2.1 


2.7 


3.6 


4.9 


6.0 


7.7 


V4 


0.95 


1.03 


1.64 


2.1 


2.5 


3.1 


5.0 


5.8 


6.3 


V5 


0.96 


1.05 


1.42 


2.4 


3.2 


3.8 


5.0 


6.7 


8.7 


V6 


0.96 


1.05 


1.44 


2.2 


3.2 


3.9 


3.5 


6.7 


8.8 


V7 


0.96 


1.05 


1.45 


2.1 


3.1 


3.9 


5.0 


6.5 


8.7 


V8 


0.96 


1.05 


1.42 


2.6 


3.0 


3.2 


5.4 


6.5 


7.4 


V9 


0.96 


1.05 


1.44 


2.2 


3.1 


3.7 


5.0 


6.3 


7.5 


VIO 


1.01 


1.14 


1.86 


2.0 


3.1 


4.2 


2.7 


5.7 


7.6 


Vll 


0.96 


1.05 


1.50 


2.7 


3.2 


4.2 


4.9 


10.5 


14.3 


Range 




0.94-1.86 






2.0-4.2 






2.7-14.3 





^ The chirp mass distribution for merging DCOs, for Zq and submodel A. The values of 
chirp mass presented are: minimum, average, and maximum in units of Mq. The range 
represents the minimum-maximum value of the chirp mass from the entire suite of models 
for each DCO type. This table corresponds to Fig. (TU] 



Table 7 

Chirp mass characteristics for Zq, submodel B^ 





min 


NS-NS 
avg 


max 


min 


BH-NS 
avg 


max 


min 


BH-BH 
avg 


max 


s 


0.96 


1.05 


1.17 


2.3 


3.2 


3.3 


5.2 


6.7 


7.4 


VI 


1.08 


1.09 


1.14 


3.2 


3.2 


3.2 


5.3 


6.5 


8.3 


V2 


0.96 


1.06 


1.53 


3.0 


3.2 


3.3 


5.6 


6.5 


7.2 


V3 


0.96 


1.05 


1.22 


2.2 


2.4 


2.6 


5.7 


5.9 


6.4 


V4 


1.03 


1.03 


1.04 




no systems 






no systems 




V5 


0.96 


1.05 


1.28 


2.4 


3.1 


3.3 


5.2 


6.7 


7.4 


V6 


0.96 


1.05 


1.43 


2.6 


3.2 


3.4 


3.5 


6.7 


8.7 


V7 


0.96 


1.05 


1.45 


2.2 


3.1 


3.5 


5.1 


6.5 


8.0 


V8 


0.96 


1.05 


1.08 


3.2 


3.2 


3.2 


5.6 


5.9 


6.1 


V9 


0.96 


1.05 


1.21 


2.2 


3.1 


3.3 


5.0 


6.2 


7.4 


VIO 


1.04 


1.13 


1.34 


2.0 


3.0 


4.2 


2.7 


4.6 


6.6 


Vll 


0.96 


1.05 


1.09 


2.9 


3.3 


4.0 


4.9 


9.1 


13.8 


Range 




0.96-1.53 






2.0-4.2 






2.7-13.8 





Same as Table [S] but for submodel B. This table corresponds to Fig. [TT1 
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Table 8 

Chirp mass characteristics for 0.1 Z©, submodel A*^ 





min 


avg 


max 


min 


r>xl-iN& 

avg 


max 


min 


OTT TI>TT 

r5rl-r>rl 

avg 


max 


s 


0.96 


1.09 


1.67 


2.1 


3.2 


6.1 


4.8 


13.2 


31.8 


VI 


1.08 


1.11 


1.56 


2.9 


3.6 


4.4 


5.9 


20.0 


32.3 


V2 


0.96 


1.09 


1.66 


2.3 


3.5 


6.5 


4.8 


17.2 


31.6 


V3 


0.96 


1.09 


1.68 


2.0 


2.9 


6.1 


4.8 


12.5 


29.7 


V4 


0.95 


1.10 


1.64 


2.1 


2.9 


5.9 


4.8 


7.6 


28.1 


V5 


0.96 


1.09 


1.66 


2.0 


3.1 


4.6 


4.8 


13.3 


31.8 


V6 


0.97 


1.09 


1.65 


1.7 


3.1 


5.2 


4.8 


13.1 


31.8 


V7 


0.96 


1.08 


1.70 


2.0 


3.1 


6.4 


4.9 


12.4 


32.0 


V8 


0.96 


1.09 


1.64 


2.0 


3.0 


6.4 


4.8 


9.0 


31.9 


V9 


0.97 


1.08 


1.66 


2.0 


3.0 


6.1 


4.9 


12.1 


31.9 


VIO 


1.10 


1.20 


2.16 


1.8 


3.4 


6.3 


2.4 


14.4 


32.0 


Vll 


0.97 


1.08 


1.61 


2.0 


3.3 


4.7 


4.8 


14.3 


41.4 


Range 




0.95-2.16 






1.7-6.4 






2.4-41.4 





Same as Table [S] but for 0.1 Z©. This table corresponds to Fig. [T^ 



Table 9 

Chirp mass characteristics for 0.1 Z©, submodel B 





min 


NS-NS 
avg 


max 


min 


BH-NS 

avg 


max 


min 


BH-BH 

avg 


max 


s 


0.96 


1.09 


1.67 


2.2 


3.1 


4.6 


4.8 


9.7 


18.3 


VI 


1.08 


1.11 


1.56 


2.9 


3.6 


4.4 


5.9 


16.1 


32.3 


V2 


0.96 


1.07 


1.66 


2.3 


3.3 


4.6 


4.8 


9.3 


18.8 


V3 


0.99 


1.12 


1.68 


2.0 


2.9 


4.5 


4.8 


6.8 


14.9 


V4 


0.98 


1.20 


1.60 


2.1 


2.7 


3.6 


5.0 


6.7 


11.5 


V5 


0.96 


1.08 


1.66 


2.2 


3.2 


4.6 


4.8 


9.8 


17.7 


V6 


0.98 


1.09 


1.65 


1.7 


3.2 


4.6 


4.8 


9.7 


18.3 


V7 


0.96 


1.08 


1.60 


2.1 


3.1 


4.6 


4.9 


9.7 


17.7 


V8 


0.96 


1.09 


1.64 


2.2 


3.0 


4.4 


4.8 


7.2 


16.3 


V9 


0.98 


1.08 


1.55 


2.1 


3.0 


4.5 


4.9 


9.3 


18.1 


VIO 


1.10 


1.20 


2.13 


2.0 


3.4 


5.0 


2.4 


10.2 


17.5 


Vll 


0.97 


1.08 


1.61 


2.3 


3.4 


4.7 


4.8 


10.6 


19.1 


Range 




0.96-2.13 






1.7-5.0 






2.4-32.3 





Same as Table [8] but for submodel A. This table corresponds to Fig. [T3l 
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Radius [R^] 



800 



Fig. 1. — The binding energy parameter, A, as a function of stellar radius, presented in lXu &: Lil ( 2010f ) (Nanjing A) for 
a typical NS progenitor in a NS-NS system (A/zams = 10 Mq) at Zq. HG stands for Hertzsprung Gap, RGB for Red 
Giant Branch, CHeB for Core Helium Burning, and AGB for Asymptotic Giant Branch. The general behaviour of A is 
described in Section! 
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Fig. 2.~ The same as Fig.[T]but for a BH-NS system (A-fzAMS = 20 M©). 
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Fig. 3. — Nanjing A of a typical progenitor for the secondary component in a BH-BH system (AfzAMS = 40 M0) at solar 
metallicity. 
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Fig. 4. — Same as Fig. [3] but for the progenitor of the primary component. 
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Fig. 5. — The distribution of progenitor (ZAMS) masses of coalescing DCOs for the standard model, submodel A, Z©. The 
top panel presents the distribution for NS-NS, the middle panel for BH-NS, and the bottom panel for BH-BH progenitors. 
Ml stands for the primary component (initially more massive, solid, blue line) and M2 for the secondary (initially less 
massive, dashed, red line). The average mass for NS-NS progenitors is 11-9 M0, for BH-NS progenitors 52-27 Mq, and for 
BH-BH progenitors 58-44 Mq (M1-M2). Note that binary evolution blurs the ZAMS mass limits for NS/BH formation 
(see Section |4J|) . 
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Fig. 6. — The distribution of remnant masses of coalescing DCOs for the standard model, submodel A, Zq. The top 
panel presents the distribution for NS-NS, the middle panel for BH-NS, and the bottom panel for BH-BH systems. Mi 
represents the primary remnant (corresponding to Mi in Fig. [5l solid, blue line) while M2 is the secondary (corresponding 
to M2 in Fig. El dashed, red line). The average mass for NS-NS systems is 1.3-1.1 Mq, for BH-NS systems 8.0-1.8 Mq, 
and for BH-BH systems 8.2-7.2 M© (M1-M2). The gap between the upper mass of NSs (2 Mq) and the lowest mass of 
BHs (5 M0) results from the use of the Rapid SN engine (see Section [2^ . 
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Fig. 7. — The distribution of chirp masses of coalescing DCOs for the standard model. The average chirp masses for 
NS-NS and BH-NS systems are ^ 1.1 Mq and 3.2 M©, respectively, for both submodels and metallicities. The average 
chirp mass for BH-BH systems, for Z©, is ~ 6.7 Mq for both submodels. For 0.1 Zq the masses are 13.2-9.7 Mq for 
submodel A and B, respectively. The maxim um chirp mass increases w ith metallicity as wind mass loss rates decrease, 
allowing for the formation of heavier BHs (see iBelczvnski et al.l ( 2010bl ) and Section ITl]) . 
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Fig. 8. — The distribution of delay times for coalescing DCOs, for the standard model. The vertical axis present the 
number of DCOs per linear time. The average delay time for all binaries is ~ 1 Gyr. 
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Fig. 9. — The distribution of mass ratios of coalescing DCOs for the standard model. The mass ratio is defined as the ratio 
of the less massive to the more massive compact object in the binary. The average values for NS-NS systems are ^ 0.85 
for both submodels and metallicities. For BH-NS binaries the average is ^ 0.22 for Zq, for both submodels and ^ 0.15 
for 0.1 Zq, also for both submodels. For BH-BH systems the average value is ^ 0.8 for both submodels and metallicities. 
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Chirp mass [Mq] Chirp mass [M^] 

Fig. 10. — The chirp mass distribution of coalescing double compact objects for all variations for submodel A and 
Z = Zq. The maximum chirp mass is found for BH-BH systems, and may reach as high as ^ 14 M©. The typical chirp 
mass for BH-NS systems is ~ 2-3 Mq , while the chirp mass for NS-NS systems peaks around ^ 1 M0 independent of the 
model. Note that the chirp masses of BH-NS systems are separated from the BH-BH values. The only exception to this 
rule is the (most-likely unphysical) VIO model, which employs the Delayed supernova engine (see Sec. 14.101 for details). 
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Fig. 11.- 
Z= Zq. 



The chirp mass distribution of coalescing double compact objects for all variations for submodel B and 
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Fig. 12. — The chirp mass distribution of coalescing double compact objects for all variations for submodel A and 
Z = 0.1 Zq. Note the dramatic increase of the maximum chirp mass with decreasing metallicity. For solar metallicity, 
the chirp mass was always below 10 Mq (Fig. [TU| . while for the majority of models shown here, the chirp mass reaches 
~ 30 M0 for sub-solar metallicity. The lack of high chirp-mass systems in model V4 is explained in See. 14.51 
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Fig. 13. — The chirp mass distribution of coalescing double compact objects for all variations for submodel B and 
Z = O.l Zq. The maximum chirp mass for submodel B typically reaches only ^ 15 Mq, as contrasted with 30 Mq for 



submodel A (see Fig. [T2|). The reason why the VI model allows for chirp mass as high as 
is explained in Sec. 



30 Mq, even for submodel B, 
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Fig. 14. — The delay time distribution of coalescing double compact objects for all variations for submodel A and 
Z = Zq. Note that especially at later times the number of sources typically falls off as cx t^^^. 
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Fig. 15.-- 



Same as Fig. [14] but for submodel B. 
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Fig. 16. — Same as Fig. [T4]but for 0.1 The average delay time for systems containing BHs decreases with decreasing 
metallicity as the components become more massive. 



35 



103 
102 
10' 
1 

10-' 
10-2 
10-3 

103 
102 
10' 
1 

10-' 

10-2 
10-3 

103 
102 
10' 
1 

^ 10-1 

^ 10-2 
\ 10-3 

T3 103 

102 

10' 

1 

10-' 

10-2 
10-3 

103 
102 
10' 
1 

10-' 

10-2 
10-3 

103 
102 
10' 
1 

10-' 

10-2 
10-3 



n I I I I [[j I I I I I M 1 1 I I I I I I M| I I I I M M| r 
















-H-Hi 








1 1 1 1 INI 


— 1 1 1 1 iiii| 



V2 



H+ 




\-HWi\ — ^ 


i' 1 1 iiii| 







+ 



V4 h 



I I l I ' M I ' MM I I N il I'll ll'lll 1 I I I M ill ^ 



V6 



ii I I I I I I 




n I I II M| I I I II I M| I I I I I I M| I I I I I I M I I 

r NSNS VI 

BHNS 

- BHBH 



il I I M il 



I 1 1 1 1 1 ; I ' 

Mill I I . Il II. 



I I I M II " I' M ll'lll I I 11 M ill 



V3 



il I I M il 1 I I M ill 1 I I M ill 1 I I llllli \ 



V5 



il I M ill 



++ 




I 

a 1 1 n il 





1 - 




1 1 1 II 1 M 


1 1 1 Mill 



11 




' I I 1_.J I I I I I I 



10' 102 104 



10» 102 103 104 

tdei/Myr 



Fig. 17. — Same as Fig. [16] but for submodel B. 
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Fig. 18. — Chirp mass — delay time diagrams for merging DCOs for submodel A, Z ~ Zq. Top panel. Results for the 
standard model. Significant populations of DCOs containing BHs are present below the 10 Gyr delay time limit (the 
horizontal dashed magenta line). Bottom panel. Results for V4. Weakly bound envelopes, represented by a very high 
A value, cause insignificant separation reduction during the CE phase. This results in only a few systems crossing the 
10 Gyr time limit and into the merging population. Both of these panels illustrate the migration of DCOs between the 
merging and non-merging populations with varying A values. The clustering of chirp masses around higher values for 
merging DCOs with BHs in the standard model comes from the fact that more massive progenitors are more likely to 
survive a strongly bound CE phase due to a larger orbital energy reservoir. Note that these diagrams do not show the 
full range of the non-merging populations; systems with delay times much larger than 10^ Myr and masses larger than 
8 Mq may also be formed. 
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Fig. 19. — Galactic merger rates from all of our models, with submodel A in the top panel and B in the bottom, for 
Z = Zq. The blue solid line re presents th e lowe r limit for predicted merger rates of NS-NS systems observed in our 
Galaxy (at 3 Myr"-'^ ) as shown in lKim et al.l ()2006() . Models yielding merger rates of NS-NS systems lower than this value 
are disfavored; these are Vl-submodel A, Vl-submodel B, V2-submodel B, V4-submodel B. Reminder: the described 
models are: V1-V4, changing A from 0.01 to 10; V5-V6, changing MNS,max from 3.0 Mq to 2.5 M0; V7 - reducing natal 
kicks for all DCOs, V8-V9, full and no natal kicks for BHs, respectively; VIO - investigating the Delayed SN engine; VI 1 
- reducing wind mass loss rates by half. 
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Fig. 20.— Same as Fig. [T9]but for Z = 0.1 Zq. 



